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heory of Stagnation Point Heat Transfer in 
Dissociated Air' 


J. A. FAY* anno F. R. RIDDELL** 


Aveo Research Laboratory 


SUMMARY SYMBOLS 


mass fraction of component 1 


The boundary-layer equations are developed in general for the 
specific heat per unit mass at constant pressure of 


ase of very high speed flight where the external flow is in a dis 


kociated state. In particular the effects of diffusion and of atom component 7, Eq. (19) 


ecombination in the boundary layer are included. It is shown €p/Cpw, Eq. (A-11 


hat at the stagnation point the equations can be reduced exactly defined by Eq. (19 
oa set of nonlinear ordinary differential equations even when the 
hemical reactions proceed so slowly that the boundary layer is 


recombination rate parameter, Eq 

hp/Epwl s 

defined by Eq. (A-14 

defined by Eqs. (A-4), and (A-7 

diffusion coefficient 

thermal diffusion coefficient 

defined by Eq. (A-13 

internal energy per unit mass of component 1 


ot in thermochemical equilibrium. 

Two methods of numerical solution of these stagnation point 
equations are presented, one for the equilibrium case and the 
ther for the nonequilibrium case. Numerical results are cor- 
elated in terms of the parameters entering the numerical formu- 


ation so as not to depend critically on the physical assumptions 

defined by Eq. (24) 

For the nonequilibrium boundary layer, both catalytic (to defined by Eq. (A-9) 
defined by Eq. (25 


recombination) and noncatalytic wall surfaces are con- 
enthalpy per unit mass of mixture, Eq. (16) 


A solution is represented which shows the transition 
perfect gas enthalpy per unit mass of component 2, 


Eq. (10) 
heat evolved in the formation of component 7 at 


the “frozen’’ boundary layer (very slow recombination 
rates) to the equilibrium boundary layer (fast recombination rates). 
A rccombination rate parameter is introduced to interpret the 


onequilibrium results, and it is shown that a scale factor is in- O°K. per unit mass 


dissociation energy per unit mass of atomic products, 
Eq. (59) 
average atomic dissociation energy times atom mass 


volved in relating the equilibrium state of a boundary layer on 
bodies of different sizes 
It is concluded that the heat transfer through the equilibrium 


stagnation point boundary layer can be computed accurately by fraction in external flow, Eq. (61 


a simple correlation formula [see Eq. (63)] and that the heat thermal conductivity 


transfer is almost unaffected by a nonequilibrium state of the 
boundary layer provided the wall is catalytic and the Lewis 


recombination rate constant, Eq. (50 
dissociation rate constant, Eq. (51 
defined by Eqs. (33), (A-5), and (A-9 
Received May 14, 1957 . = ee ee 
t This work was sponsored by the Ballistic Missile Division, E Di" pte/%, eevee Loews Newber 
ARDC, USAF, under Contract AF 04(645)-18. Many mem ‘ Lewis Number for atom-molecule mixture 
bers of the staff of the Avco Research Laboratory have con- 
tributed to this paper. In particular, thanks are due to A. R 
Kantrowitz, who initiated the work, and to N. H. Kemp, 
who meticulously reviewed the results and suggested several 


Number near unity. 


molecular weight of molecules in atom-molecule 


mixture 
number density, moles per unit volume 
Nusselt Number, Eq. (43) 


Improvements. pressure 


In carrying out the computations the work of D. Goldberg, = heat flux 
R. D. Laubner, J. Levin, and W. M. Wolf is gratefully acknowl- ; ’ . 
s : vector mass velocity 
edged 
* Consultant, Aveo Research Laboratory, and Professor, De- ‘ = vector diffusion velocity, Eq. (2 


partment of Mechanical Engineering, Massachusetts Institute of cylindrical radius of body 


body nose radius 


Technology A 
= Reynolds Numbers, Eq. (44 


** Principal Research Scientist. 
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R,; = gas constant of component 7 

Re = gas constant of mixture 

R = universal gas constant 

s defined by Eq. (27) 

7 = absolute temperature 

Ty = vibrational temperature for molecules, taken as 800° 
K. 

u = x component of velocity 

v = ycomponent of velocity 

W; = mass rate of formation of component 7 per unit 
volume and time 

x = distance along meridian profile 

y = distance normal to the surface 

ay, a defined by Eq. (A-5) 

Bi, Bo defined by Eq. (A-7) 

v1, ¥2 defined by Eq. (A 6) 

n defined by Eq. (22) 

0 defined by Eq. (26) 

Ib = absolute viscosity 

v = kinematic viscosity 

é defined by Eq. (23) 

p = mass density 

o = Prandtl Number ¢,u/k, taken as 0.71 in numerical 
calculations 

? = dissipation function 


Subscripts 


A = atom 

1 = ith component of mixture 

e = external flow conditions 

uw ‘= wall 

s = stagnation point in external flow 
M = molecule 

m = mixture 

E = equilibrium 


(1) INTRODUCTION 


< PROCESS of aerodynamic heat transfer at hyper- 
sonic velocities is complicated by two features not 
normally present at low velocities. The first is the 
possible dissociation and ionization of air due to high 
static temperatures encountered where the air is de- 
celerated by shock waves, by viscous forces in the 
boundary layer, or at a stagnation point. Because dis- 
sociation and ionization (and their reverse processes, 
recombination) proceed at finite rates, thermochemical 
equilibrium is not necessarily achieved throughout the 
flow field, and such rates are therefore an essential in- 
gredient of the flow process. Secondly, diffusion of 
atoms and ions, which subsequently recombine with a 
high specific energy release, may appreciably add to the 
heat transferred normal molecular conduction. 
While there may be other physical phenomena also 


by 


present, such as radiative effects, it is the purpose of the 
present analysis to include only the two effects pre- 
viously noted in an otherwise classical viscous flow 


problem. 
Some aspects of this problem have already received 
attention. Moore! considered a dissociated laminar 


boundary layer on a flate plate in air with a local com- 
the thermochemical eaui- 
1.e., a recombination rate constant sufficiently 
Hansen* noted 


position determined by 
librium 
great to maintain local equilibrium. 
that Moore had miscalculated the Prandtl Number for 


dissociated air and an analysis similar to Moore's but 


‘ 
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following Hansen’s suggestion was made by Romig “a 
Dore.* Beckwith‘ considered the heat transfer to the 
stagnation region of a blunt nosed body, using integral | 
methods.* Crown® the stagnation 
point problem, using a modified Crocco method for 
solving the boundary-layer equations. Finally, Mark‘ 
treated the stagnation point equilibrium boundary 


also considered 


layer with variable fluid properties.T 

In all these analyses no detailed distinction is made 
between the roles of atomic diffusion and molecular 
conduction in transporting energy to the wall.f If one 
considers the energy transport through a motionless 
dissociated gas with temperature and concentration 


gradients, the energy flux is approximately 
g = k grad T + ha°Dp grad cz, 


where k& is the ordinary thermal conductivity, 7 is the 
temperature, 44° is the dissociation energy per unit 
mass of atomic products, D is the atomic diffusion co- 
efficient, p the density, and c4 the atomic mass fraction. 
The first term is the usual transport of kinetic, vibra- 
tional and rotational energy and the second is the 
transport of potential (recombination) energy. Even 
in the “equilibrium’’ boundary layer, the latter term 
should be taken into account in the energy equation, 
and will constitute a significant contribution where the 
atom concentration and diffusional velocities are 
noticeable. 

These distinctions were first pointed out by Fay." 
Subsequently, Lees’ considered in detail the laminar 
heat-transfer problem in dissociated air, including the 
effects of atomic diffusion, and suggested several ap- 
proximations to facilitate the solution of the boundary- 
layer equations. He considered the limiting extremes 
of the recombination rate constant which (as discussed 
below) lead to simpler solutions than the general case, 
and suggested expressions for the heat transfer for both 
cases. Similar arguments, but in less detail, are ad- 
vanced by Kuo.'* 

Apart from the physical mechanism of heat transfer, 
there are several relevant aerodynamic considerations, 
the foremost of which is the question of shock-wave 
boundary-layer interaction at a sharp leading edge. 
From a mechanical point of view, it does not appear 
possible to maintain sharp leading edges with their 
attendant high heat-transfer rates, so that a finite 
radius of curvature appears mandatory. Under these 
conditions a distinct layer will exist independent of the 
detached bow shock wave if the boundary-layer thick- 
ness is much less than the shock detachment distance. 


* It is noted that the variation of viscosity through the bound- 
ary layer, which materially affects the heat-transfer rate, cannot 
be accounted for by the integral method. 

+ For a more detailed discussion of numerical results, see Fay, 
Riddell, and Kemp. !* 

t Crown,5 for example, assumes that the potential energy of 
dissociation (which is actually carried by diffusion of atoms) 1s 
transported in the same manner as the internal energy of the 
molecules, which implies equal diffusivities of molecules and 
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THEORY OF 


Since the former varies as the inverse square root of the 
Reynolds Number and the latter is independent of 
} Reynolds Number, there is a minimum Reynolds Num- 
? ber below which the shock wave and boundary layer 
However, since low Reynolds Numbers at hy- 


—~—— 


merge. 
personic velocities can only be attained by reducing the 
density (and hence increasing the mean free path), this 
limit may not be reached before free molecule flow en- 
sues. It thus appears that stagnation point boundary 
layers are cuite relevant to hypersonic continuum 
flows. 

Boundary layers at locations other than the stagna- 
tion point are also of interest, but will not be considered 
herein in detail. There are certain difficulties inherent 
in such solutions when the recombination rate is finite 
which will be discussed more fully below. 

For any boundary-layer calculation, the behavior of 
the free stream outside the boundary layer must be 
known. If the boundary-layer thickness is small com- 
pared with the nose (and also detached shock wave) 
radius of curvature, then only fluid which has passed 
through the normal shock wave close to the axis (or 
plane) of symmetry will enter the boundary layer, and 
this constitutes the ‘‘free stream.”’ If it is assumed that 
this stream has attained thermochemical equilibrium 
by the time it reaches the stagnation point, then the 
iree-stream conditions for the stagnation point bound- 
ary layer are those of an equilibrium gas. This would 
appear to be the case most likely to be encountered.* 

Whether the free stream in regions other than the 
stagnation point remains at equilibrium again depends 
upon the recombination rate since, in flowing to regions 
of lower pressure, the gas is expanded and cooled. 
Again, there will be no general solution which will in- 
clude all possible situations. 

The two extremes of recombination rate give rise to 
simpler solutions than the general case. For suf- 
ficiently small recombination rate, the concentration of 


| atoms (or ion pairs) is determined by the diffusive flow 


from the free stream to the wall where recombination 
would occur, and would bear no relation to the thermo- 
chemical equilibrium concentration corresponding to 


the local temperature. In such a ‘‘frozen’’ boundary 


layer, the temperature and concentration distributions 


are practically independent of one another. On the 
other hand, for a sufficiently large recombination rate, 
constant thermochemical equilibrium would prevail 
throughout and either the temperature or concentration 
distribution is a sufficient description of the thermody- 
namic state of the boundary layer. It is, therefore, to 


be expected that the distribution of atoms in the 


*The time for the gas to come to equilibrium behind the 
normal shock wave depends upon the kinetics of the dissociation 
process. It is merely noted here that immediately behind the shock 
wave, before dissociation begins, the translational temperature is 
extremely high and thereby promotes high dissociation rates 
Even at high temperatures equilibrium may not be attained soon 
enough if the density is sufficiently low, but this would probably 
oecur only near the free molecule flow regime. 
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“frozen”? and boundary layers will be 
quite different 


Despite the importance of the process of dissocia 


“equilibrium” 


tion and recombination in determining the thermody 
namic state of the air throughout the flow field, their 
effect on heat transfer is secondary. This is most easily 
seen by again considering the heat flux through a stag 
nant gas as given above, and replacing grad 7° by 


(1/c,) grad h, where h is the perfect gas enthalpy. If we 
also make the approximation that Dpc,’k (Lewis 


number) is unity, the heat flux becomes 
g = (k/c,) grad (h + caha°) 


that is, the heat flux is determined by the chemical en- 
thalpy (perfect gas enthalpy plus enthalpy of forma 
tion) difference between free stream and wall. Whether 
atoms recombine in the boundary layer or on the wall 
makes no great difference since the energy is conducted 
about as readily by normal conduction as by diffusion 
when the Lewis Number is approximately one. On 
the other hand, if the heat transfer is to be known more 
exactly, then such relevant effects as variation of trans- 
port coefficients with temperature, variation in heat 
capacities, actual Lewis Number, etc., must be ac- 
counted for properly. It will be seen that it is the in- 
fluence of these effects which constitutes the principal 
departure from an extrapolation of the classical theory. 

Since a considerable fraction of the heat may be 
transported by atomic diffusion toward the wall fol 
lowed by recombination on the surface, it would be 
possible to eliminate this fraction of the heat transfer by 
using a noncatalytic surface. However, such a scheme 
is useful only if the atoms do not first recombine in the 
gas before reaching the wall. The flight conditions 
under with the gas phase recombination is slow enough 
to permit atoms to reach the wall, and the resulting heat 
transfer with both catalytic and noncatalytic surface, 
have been determined for the stagnation point flow. 

For a discussion of experimental techniques and re- 
sults of shock tube measurements of stagnation point 
heat transfer in dissociated air, see reference 15. 


BOUNDARY-LAYER EQUATIONS IN A 


DISSOCIATED GAS 


(2) LAMINAR 


The general equation of continuity for any species 7 is 


div} p( g + g ie it wW (1) 
where p is the mixture density, c; is the mass fraction of 
species 7, w,; is the mass rate of formation of species 7 
per unit volume, g is the mass averaged velocity and q; 
the diffusional velocity of species 7 measured with re- 
spect to g. gi may be given by 


> 


Gi = grad 7° (2) 


/ ry 
—(D c;) grad c; — (D iS 
where D; and D,' are, respectively, the molecular and 
thermal diffusion coefficients of species 7, and 7 is the 


temperature. The first term on the right of Eq. (2) is 
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due to concentration diffusion and the second is due to 
thermal diffusion (pressure diffusion is neglected). 

The use ofthe gradient of mass fraction instead of 
mole fraction is particularly useful since, in the case of a 
bimolecular mixture, D; becomes the bimolecular diffu- 
sion coefficient ();.) which is practically independent 
of composition. For multicomponent mixtures, the 
diffusional velocity of any one component depends, in 
general, upon the concentration gradients of all the com- 
ponents, and a method for determining the diffusional 
velocities has been suggested by Hirschfelder, Curtiss, 
and Bird.* In such a case Eq. (2) is not exact (except 
for equal diffusivities and molecular weights of all com- 
ponents), but is a useful approximation. <A dissociated 
gas (such as air) in which all the molecules (or atoms) 
have nearly the same molecular weight and probably 
similar collision diameters may be considered primarily 
a two-component mixture with atoms and molecules as 
species, for which Eq. (2) is adequate. 

In order to conserve mass in any chemical change it is 
necessary that 


2w, = 0 (3) 


and from the definition of the mass averaged velocity it 


follows that 
Lai, = 0 (4) 


Hence, Eq. (1) summed over all species gives the usual 
form for the continuity equation, 


div (pg) = 0 (5) 


The energy equation may be written for a moving 
element of fluid 


pq -grad(Xc,e;) = div(k grad 7) — 
div(Zp q icih;) + Ywh? + p div q + (6) 


where ¢;, h;, and h;® are, respectively, the specific in- 
ternal energy, enthalpy, and heat of formation of 
species 7, k is the thermal conductivity for transport 
of kinetic, rotational, and vibrational energy, p is the 
pressure, and ® is the dissipation function. 

The first term on the right of Eq. (6) is the internal 
energy increase due to normal heat conduction, the 
second that due to fluid diffusing across the boundary 
of the element, the third that due to chemical reaction, * 
and the fourth and fifth that due to the work of the 
pressure and the viscous forces, respectively. 

It will be assumed that the gas is a mixture of perfect 
gases, so that for each component 


Pi = pR,T (7) 
where RX; is the gas constant for species 7, and for the 
mixture 

pb = eR, T (S) 
where the gas constant for the mixture (R,,) is given by 


* The heat of formation h;° may be taken as zero for the 
molecules and negative for the atoms, in which case the enthalpy 
h {see Eq. (16)] is always positive. 
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Rh = Zea; (9 ] 


It then follows that the enthalpy and internal energ 


are related by 
h,=e,+ RT If 


By combining Eqs. (5) and (7) through Eq. (10) wit 
Eq. (6), the steady-state energy equation reduces t 


pq-grad(2c;h;) = div(k grad T — 
Lp qicihi) Lwh? + q ‘grad S+- 2 (i 


A further simplification occurs by combining Eqs. | 


and (11) to eliminate the term Yw,h ;’: 


pq-grad| Sc(h; — h)} = div{k grad T — 
Legiedhi —h?)} + g ‘grad pt+e® (12 


By making the usual boundary-layer assumptions 
Eqs. (1), (5), and (12) reduce to the following form for 
body of revolution, if the boundary-layer thickness j 
small compared with the radius of curvature and centrii 


ugal forces are neglected :T 
(pru), + (prv), = O (13 
pucir + pvcy = | Dipciy, + Di" pciT, Ti, +w; (14} 


puh, + puh, = (RT,)y + upe + 
u(uy)? + {| LDip(h; — hci, + 
DD" pei(h; — h®)T, T}, (15 


where ¢ is the radial distance of the body surface fron 
the axis of revolution, u and v are the velocity com 
ponents in the x and y directions (tangential and normal] 
to the surface, respectively), u is the absolute viscosit) | 
and the enthalpy for the mixture is 





h = 3c,(h, — h;) (16 
The corresponding equation of motion is 
puu, + puuy = —pr t+ (uuy,); (17! 


Since the transport coefficients are, in general, fanel 
tions of temperature and composition, it may be de} 
sirable to use 7 rather than / as the dependent variabk; 
in the energy equation. Noting that h; is a function 0} 
temperature alone, we have 
grad h =| c,(dh,;/dT)} grad T + | 

Dh; — hi®) grad c; (ist 


Letting 
Cp = 2c,(dh,/dT) = Zep: (1 
where c,; is the constant pressure specific heat fo 


translation, rotation, and vibration, the energy equa 
tion (15), in combination with Eq. (14), becomes 


| pul + poT,} = 
(RT,), + ups + u(u,)? rwi(h® — hy) + 
Lop Diplyy + Dit pc;T,/T)T, (2 
For a gas at equilibrium, on the other hand, it is foun 


{ Subscripts x and y (or £ and 7 below) are used to denote pa’ 
tial differentiation. 
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| more convenient to use the enthalpy as a variable rather 


than temperature, in which case by combining Eqs. (17) 


and (1S) with Eq. (15), we obtain 
ou(h + uw? 2), + pr(h + u? 2), = \ (Rk, e,) X 
h+ uw?/2),5, + (C 2)[a — (& @,) (u?),5, + 
i=(Dip — (k E,) (hy — AY )c, 4 
Z(DT pe:/T)(h, — h)T,f, (21 
Eys. (13), (14), (17), and (20) or (21) thus constitute 


the system whose solution is required. 


3) SIMILAR SOLUTIONS 


As is usual in boundary-layer problems, one first 
seeks solutions of restricted form which permit reducing 
exactly the partial differential equations to ordinary 
differential form. An easily recognizable case is that of 
the stagnation point flow, where, because of symmetry, 
all the dependent variables are chosen to be functions of 
yalone, except « which must be taken proportional to x 
times a function of y. This also appears to be the on/y 
case for which exact ordinary differential equations 
may be obtained regardless of the recombination rate. 
For the flat plate and cone, exact solutions exist only 
for the extreme values of the recombination rate con- 
stant, that is 


“frozen”’ or in the thermodynamic equilibrium. 


when the boundary laver is either 

For all other cases, certain degrees of approximation 
are required. Following standard procedures, one first 
tries solutions for which the velocity and enthalpy pro- 
files remain similar to themselves, at least for an ap- 
the Such “‘locally 
similar’’ approximate solutions may be obtained for the 


preciable distance along body. 


‘frozen’’ or ‘‘equilibrium’’ boundary layer, and also for 
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The existence of a finite recombination rate is not the 
only hindrance to obtaining exact solutions in regions 
other than the stagnation point. In general, the varia- 
tion of the velocity, thermodynamic variables and trans- 
port coefficients in the free stream and along the wall of 
the body of arbitrary shape preclude exact solutions 
(except for the cone and flat plate). An important 
aspect of the locally similar solutions is a proper ac- 
counting for these variations such that the heat transfer 
may be determined for any point of the body of arbi- 
trary (but regular) shape. Although this problem will 
not be treated in this paper, we will start with a trans- 
formation suitable for locally similar solutions. 

With this in mind, we choose the following transfor- 
mation of the independent variables « and y, which in 
cludes the usual Howarth and Mangler transformations, 
as proposed by Lees.’ 


pdy (22)” 

JZ j 

= } Probhiclt FUN (23)" 
{ 


where u, is the velocity at the outer edge of the bound- 
ary layer. In addition, the following dimensionless de- 


pendent variables are chosen 


Of/0n = u/u.; f = | (Of On)dn (24 
JV 

g= (hk + u?/2)/h (25 

6= T/T, (26 

‘ = € ( (27 


where the subscript e refers to values in the local free 
stream, and s to values in the free stream at the stagna- 


an arbitary recombination rate in certain restricted ; : i oo dena papde: i 
in tion point. Substituting in Eqs. (13), (14), (17), (20), 
cases. 
and (21), there results 
pu = — (VY 2E fe + I/ VY 2Er + VY 2E fynz] (28) 
r\(Lisip +L ’s 6, 0)|, + fsi, + [REw; pu.c;(dé dx)| = oe si — feStq) + 2f,5 [d(In « d(In &)] (29 
(fon)» + ffom + 2[d(n u.)/d(In €)][(p. p) — f,7] = 2S afr — Sehr (30) 
(@, /é a)6.|, + (Co/Cen lite > [2Ew pu,(d& dx) |[(h o— h,;) ow + (u,? Cone 10 ee + 
DhGos/ Ci tt o)(L,S;, + L,'s6,/0)0, = f,, 2(Ep Cpe)O|d(In 7.) d(in &)] 4 
(21.7 Enel) (pep) [d(n u,) ‘d(In E)]} + 2G, Ene d& 0: fO, (31) 
((l/o)e,], + fz, + (ue hy — a-Yif, fort, + 
IU a)S(ce hy (hy — AY) (Li, — 1)s;, + L,7s 6,6]; , RE(f ee — Seg, 3) 
where i= PL Pwby 


and where o and L; are the Prandtl Number (¢,u &) and 
Lewis Number (Dé, k), respectively. The subscript w 
refers to values of the variables at the wall. 

For the dependent variables f, g, 6, and s, to be func- 
tions of 7 alone, it is first necessary that the thermody- 
namic state variables be unchanging in the free stream 


and along the wall as € increases, a condition satisfied 
at a stagnation point, or along a cone or flat plate. It 


is further necessary that the source term—.e., the term 


involving w,—in Eqs. (29) and (31) depend upon 


alone. Since the chemistry of the recombination proc- 
ess is believed to depend only upon the local thermody- 
and 


namic variables—that is, w; is a function of p, 7, 


s;—this second condition may be satisfied in three dif- 
ferent ways: (a) w; = O (frozen boundary layer), (b) 
u,d(In £) dx = constant, which is satisfied at a stagna- 


tion point, or (c) the recombination rate is sufficiently 


* For two-dimensional flows let r = constant 
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large to maintain thermodynamic equilibrium,* in 


which case w; is determined through Eq. (29) from the 
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solution of Eqs. (28), (30), and (32). 


rs) 


For the stagnation point, then, Eqs. (29) to (32) r 


duce to 
1 ae Fame ‘ . P : , 
[(l/o)(Lis;, + Li 5:T,/T)], FSi T ) 2(du,/dx),¢ —'(w;/ pcis) = O 34 | 
y 9\ 5 ol ian 
(fan) n TOSS (1/2)\ (psp) — Jeni = 9 29) 
(€,/Epw)(l/ a) |, t+ (Ep/Cpw)fO, + ) 2(du,/dx),¢ 12(w,/p)[(hY — hj) Cpwl 3) + 
D(Cyi/Cow) (Cigl/o)(LiSig + Ly 58,/0),0, = 90 (3 
\(/o)g, J, + fg, + ((l/o)Z[ci(hi — h®)/h,] \(L; — 1)s;,, + L's, 0), = 0 37 
where it has been assumed that u* < h,. 
For a ‘locally similar”’ solution away from the stagna- or 
tion point it is required that the terms on the right-hand 7 
gq = |(k €,)(Oh Oy) ]y-0 + [2(k C,)(hi — hi’) X 


side of Eqs. (29) through (32) all be negligible compared 
with those on the left, and thus the equations may be 
integrated with respect to 7. The dependence upon & 
of f, g, #, and s; is thereby implicitly determined by the 
variation of p, 7, s;, u, and the transport coefficients in 
the free stream and along the wall, this variation being 
determined by the aerodynamics of the external flow. 
For Newtonian flow the distribution around the 
body of the relevant quantities is easily found, all quan- 
tities except du, dx being determined by the local body 
slope— i.e., local pressure. Since du, dx depends upon 
the pressure gradient and hence body curvature in the 
meridian plane, there is introduced an extra parameter 
in the locally similar solutions. Lees’ points out that in 
the equation of motion (30) the pressure gradient affects 
only a single term, that this term is generally small, and 
that neglecting this term in the equation of motion will 
This 


appears to be a greatly desirable simplification in de 


only slightly affect the solutions for g, @, and s;. 


termining the heat transfer at other than the stagnation 
point. 

The use of locally similar solutions can only be justi- 
fied a posteriori in each particular case by determining 
from such a solution the magnitude of the terms in 
Eqs. (29) through (52) which were neglected. Pre- 
sumably an iterative scheme could be devised to im- 
prove on such solutions. It does not appear possible to 
formulate a general criterion for determining the limit 
of applicability of the locally similar solutions except a 
loose physical argument that conditions in the free 
stream and along the body must change only slightly in 
a distance of many boundary-layer thicknesses. 


(4) HEAT TRANSFER 


The local heat-transfer rate to the body g is deter- 
mined by the sum of the conductive and diffusive trans- 
ports, the latter being included only when the atoms 


recombine on the wall. Thus, 


— h;®) x 


7 = [R(OT Oy) |, 0 + [Sp(h 
T)(oT oy); l, 0 


|D(Oc; Oy) + (D,"e (38) 


* Strictly speaking, the fluid cannot be exactly in thermody 
namic equilibrium, for if this were so, there could be no net rate 
of change of composition following a particle path. For large 


enough rate constants, the fluid will be very close to equilibrium. 





}(L; — 1)(0c¢; dy) + (L;"¢;'T)(OT /dy)} Jy (39 


In terms of the dimensionless temperature and en 
thalpy distributions, this becomes 


Gg = (rkuputtel e/V 2£)}0, + Deir X 
[(h; — hy) /é,T.\(Lisin + L;'s0,/0)),-0 (40 


or 
Gg = (Repth V/ 2E Enw) X 
1g, + Le,[(hi — hi®) h,] X 
((L; — 1)s;, + L,'s6,/0]},-0 (41 
For stagnation point heat transfer we note that 
rpytt. YW 2E = \(2/v,)(du, ‘dx) baie $2 


where y, is the kinematic viscosity at the wall. 


It is also possible to define local Nusselt and Reynolds 


Numbers based on the local coordinate x. Thus, de 
finingtT 

Nu = GxCpw/Ruylhs — h 13 

Re = u,X/V 14 


the heat transfer at the stagnation point may be written 
as 
g = (Nu V Re)V Puktuc(due dx). (hy — hy) /o] (45 


where 


6, + Zcis((hi — hi)/CpT;] X 
(Lisip + Li’ S8,/0)t,-0 (46a)t 
Nu/~/ Re = [VW 2/(1 — gw)] X 
12, + Deis((hi — h3°)/hs| X 
((L; — 1)s;, + L,'s6,/0]},-0 (46d)! 


[t should be emphasized that the heat transfer is in- 
dependent of the particular choice of reference density 
and viscosity (p,,“,) Which appears in the definition of 


Eq. (23). The use of wall values of p and yu for this pur- 


+ It isclear from the expressions (40) and (41) that the enthalpy 
difference between wall and free stream, not the temperature 
difference, is the proper ‘‘driving force’’ for heat transfer 

t For conical or flat plate flows, the factor 1/2 in Eq. (46 


9 


\/3/2 or 1/+/2, 
reference dimensions for Nu and Re are the distances from the 


should be replaced by respectively, when the 


apex or leading edge, respectively, measured along the surface 
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pose is merely a convenience, but a logical choice in that 
| [Eq. (33)] at the wall (y = 0). 
However, another choice, such as free-stream p and uw, or 


it ensures that / = 


evel pu equal to an arbitrary function of xv, could also be 
used. In each such case, however, the solutions of the 
boundary-layer equations give different g(n), @(m), ete. 
It is clear from an examination of these equations that 
each such solution may be transformed into another by 


a simple change of scale of n and f defined by 


V Pwo 2 = V Poko 1 (47) 
V Pwltw f = V pomo fi (48) 

and 
g(n) = go(V/ poo Puktw 10), etc. (49) 


and where pouy is the new reference value used in defining 
m. AS a consequence the heat transfer [see Eqs. (38) 
and (39)] can be seen to be independent of the choice of 
Pole 

On the other hand, if one attempts to estimate g,(Q) 
for the case of pu not constant from a solution for pu 
constant,* which is the problem considered by Lees,’ 
the reference state (poo) is quite relevant since by its 
proper (and arbitrary) selection the wall gradient of g 
may be made to equal numerically that for the case of 
pu constant. A discussion of this point has been given 
by Probstein,? and by Fay, Riddell, and Kemp.'” 


(5) TRANSPORT PROPERTIES 


To obtain a numerical solution of the boundary-layer 
equations it is necessary to know the viscosity, Prandtl 
Number, Lewis Number, and thermal diffusivity as 
functions of the local temperature, density, and com- 
position. For the high temperatures under considera- 
tion, these properties have been estimated,f the principal 
lack of knowledge of 


atom-molecule interaction potentials. 


uncertainty resulting from 
The particular 
results of this study which are pertinent to the bound- 
ary-layer calculation may be summarized as follows. 

a) The viscosity of equilibrium air, determined by 
assuming that all particles possess the same interaction 
potentials, does not vary more than 10 per cent from 
Sutherland’s formula below 9,000°K. 

b) Both Prandtl and Lewis Numbers do not change 
appreciably with temperature (below 9,000°K.), the 
value of the latter being uncertain, but estimated as 
about 1.4. 

In light of the uncertainties in viscosity and Lewis 
Number, it is clear that the calculated heat transfer 
using these estimates is also uncertain. However, the 
numerical solutions of the boundary-layer equations 
may be obtained for a range of possible Lewis Number 


and viscosity variation, such solutions being generally 


valid when expressed explicitly in terms of the property 


\s, for example, the solutions of Cohen and Reshotko.! 
t These estimates were made by Dr. S. Penner of California 
Institute of Technology and M. M. Litvak of Cornel! University. 
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variations. Thus for the numerical solutions, the 
Prandtl Number was held fixed at 0.71, the Lewis Num 
ber was assumed constant through the boundary layer 
at values ranging from one to two, and the viscosity 
variation was determined by Sutherland’s law for the 
equilibrium air and as calculated by Penner and Litvak 
for nonequilibrium composition. { 

It is easily shown that thermal diffusion is unimpor- 
tant for the equilibrium boundary layer at stagnation 
than 10,000°K. For the frozen 


boundary layer this is no longer necessarily the case; 


temperatures less 


however, thermal diffusion was neglected in all the 


numerical solutions reported in this paper. 


(6) RECOMBINATION RATE 


As suggested by Davidson,** the recombination rate 
of atoms is determined by a three-body collison: 
atom + atom + particle ~ molecule + particle 
For this process the rate of disappearance of atoms may 
be written as 
dN. ‘dt —K,Nx?NT™- 50 


where the temperature dependence is as suggested by 
to be 5 X& 10 
' for oxygen recombination when 7 


Davidson, who also estimated A,7 
ce.” mole~* sec. 
300°K.77 N4 and are the number of moles per cm. 
of atoms and particles, respectively. 

Atoms will be produced by the reverse of the above 
reaction, so that the net rate of production may be 
written as 


dN, dt = —K,Na’NT- 6 + KioT)NuN (51 
Since under equilibrium conditions, there is no net pro- 
duction of atoms, A.(7) may be solved for in Eq. (51 
in terms of the equilibrium atom and molecule concen 
trations which would exist at the local temperature 7 
and pressure. Substituting this value of A» in Eq. (51 


there results 


dN4/dt —K,N4’NT~'* X 


; l \ Na Nw} Na 


where the subscript £ refers to the concentrations which 
would exist at thermodynamic equilibrium at the local 
pressure and temperature, and Vy is the number of 
moles of molecules per em 

Considering a gas mixture of diatomic molecules of 


t Results of the experiments of reference 15 indicate good 
agreement with the present estimates of transport properties with 
Prandtl Number 0.71 and Lewis Number i.4 

** Dr. Norman Davidson, California Institute of Technology 
private communication 

tt Recent experiments reported by J. Camm and J. Keck (see 
Bulletin of the American Physical Society, Series II, Vol. 2, No. 4, 
p. 216, 1957) indicate that the recombination rate at 6,000°K 
may be higher than that obtained by extrapolating Davidson’s 
estimate to such a temperature. There is nothing known con 
cerning this rate at the low wall temperatures pertinent to the 


boundary-layer problem 
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molecular weight .J/ and atoms of weight .\/ ‘2, the atom 
mass fraction and total density are, respectively, 


Ca = Na (Na + 2Ny) (53) 
and = M(Nuy + Na 2) (54) 


The net mass rate of formation of atoms (w.) may now 
be determined in terms of these quantities to give 


ws = (M/2)(dN4,/dt) | 
— (2K,p*T—'!-*/M?) X (59) 

[((1 + c4)(€4? — Cag?)/U — can?) |) 
Neglecting the term (1 — c¢4,"), which does not vary 


by 25 per cent in the cases considered, the relevant 
source term in the diffusion equation (34) for atoms re- 


duces to 


1 2(du, dx), —(w;'p) = —[Kip2T,-**/R2(du, dx),|] X 
10-3:5(¢42 — can?) (1 + 4)} (56) 


where ® is the universal gas constant. The coefficient 
of the term in braces on the right of Eq. (56) contains 
all the dimensional effects of the finite recombination 
rate problem, and will be termed the recombination rate 
parameter (C\)—.e., 


I (57) 


C, = Kip/T,-**R-*(du,/dx), 
The term (du, dx),~' is approximately the time for a 
particle in the free stream to move a distance equal to 
the nose radius, and thus also the time for a particle to 
diffuse through the boundary layer at the stagnation 
point. The factor which multiplies it in Eq. (57) is the 
reciprocal of the lifetime of an atom, so that the recom- 
bination rate parameter is the ratio of the diffusion 
time to the lifetime of an atom. Because the diffusion 
time contains the body nose radius while the lifetime 
does not, a scale effect is introduced by the chemical 
change which is not accounted for in the Reynolds 
Number. Thus, similar flows require equal Reynolds 
Numbers and recombination rate parameters. 

A term similar to the left side of Eq. (56) appears in 
the energy equation (36) except that it is multiplied by 
a dimensionless dissociation energy. A part of this 
term, Yw,h;, is negligible since 2w; = O and h; per unit 
mass for vibrationally excited molecules and atoms are 
nearly equal, being in the ratio 9/10. 


(7) NUMERICAL SOLUTIONS 


The boundary-layer equations derived above [Eqs. 
(34) through (37)] can be put in a form suitable for 
numerical integration. The details of how this was 
done are given in the Appendix. Numerical solutions 
were obtained by use of an IBM 650 digital computer. 

Many combinations of the parameters o, L, and C, 
are possible, as well as extreme ranges in the free-stream 
(stagnation point) and wall boundary values. It was 
therefore decided to restrict the solutions to the follow- 
ing values: 


(a) «¢ = 0.71 
(b) L = 1.0, 1.4, 2.0; L’ = 0 
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(ec) Oo & Ly. Ss 
d) Stagnation point conditions corresponding to ther 
modynamic equilibrium at velocities between 5,S0( | 
ft. per sec. and 22,800 ft. per sec. and at altitudes oj 
25,000 ft. to 120,000 ft. (Equilibrium air properties 
were determined from the tables prepared by th 
National Bureau of Standards 
e) Wall temperatures from 300° k. to 3,000°K. 
Wherever possible, only one parameter or boundary | 
condition was varied in a series of calculations in order 
to determine its individual effect. It was not believe 
necessary to establish such effects for all possible com 
binations of the remaining parameters. 
the values of the parameters and the boundary condi 


tions for each computation may be obtained directly | 


from the authors. 

The heat-transfer parameter (.Vu vy Re 
from each integration depends not only upon. the 
the dependent 


determined 


parameters and boundary values of 
variables involved, but also upon the assumed variatiot 
of viscosity with temperature and composition. In ad- 
dition, further simplifications as explained below for the 
individual cases were also made. In order to make the 
results less dependent upon the specific assumptions 
made, the heat-transfer parameter was numerically 
correlated with the parameters and boundary values in 
what seemed to be a suitable manner. It is believed that 
not too widely different assumptions would give results 
which would fall within the same correlation. 
S) THE EQUILIBRIUM BOUNDARY LAYER 

A first series of equilibrium boundary lavers was 
computed by solving Eqs. (35) and (37) simultaneously 
(hereafter denoted as Method 1). The solution requires 
specifying p, p and / = py py, as functions of g 
typical variations of these quantities being shown in 
Fig. 1. The calculated values at various altitudes of 
interest are compared with the fitted curves used in the 
computing program, from which it was concluded that 
there was a negligible effect of altitude variation on these 
Specific details of the method of calculation 
It was found that moderate 


functions. 
are given in the Appendix. 
changes in the distribution for identical end values re- 
sulted in negligible changes in the heat-transfer par- 


ameter. 
For L; = 1, the equations are similar in form to those 
solved by Cohen and Reshotko,'® and become identical | 


at low enough stagnation temperatures when 
proximately constant. Solutions were obtained for the 
range of velocities and wall temperatures given above 
and the heat-transfer parameter was found to depend 
only upon the total variation in py across the boundary 
layer, in accordance with the relation 
Nu/vV/ Re = 0.67 (petts/ Putte) ®*4 (58 
The numerical correlation leading to Eq. (5S) 1s shown 
in Fig. 2. The solution of Cohen and Reshotko!® for 
/ = 1 is also plotted, after correcting for Prandtl Num- 


0.4 


ber by multiplying their result by (0.71 


Tables giving 
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Fic. 1 Variation of pu/ pu. and p, p through a stagnation 


point boundary layer, shown as a function of enthalpy ratio g. 
Stagnation conditions are for a flight velocity of 19,700 ft. per 
sec. at various altitudes. 7, = 300°K 


An alternative procedure for the equilibrium bound- 
ary layer is to solve Eqs. (34), (35) and (36) simultane- 
ously (hereafter denoted as Method 2). An appreciable 
simplification results if we consider air to be composed 
“air’’ molecules and “‘‘air’’ 


of only atoms having an 


average heat of formation given by 


sa, = > c;.(—h;) pe (59) 
atoms atoms 


where the summation extends over atomic oxygen and 
nitrogen only. Thus only one diffusion equation (34) is 
needed for the diffusion of air atoms. With this simpli- 
fication the equilibrium boundary layer may be treated 
by eliminating the term involving w,; between Eq. (34) 
and (36). A solution is then possible when ¢, and / are 
specified as functions of s4 and 6, and s4 is specified as 
a function of @ through the known equilibrium atom 
(Details of the 
approximations made are given in the Appendix.) 
This alternative solution was found to give very 


fraction as a function of temperature. 


closely the same results as the Method | for a Lewis 
Number of unity, and the results are compared with 
Eq. (5S) in Fig. 2. For a Lewis Number of unity, 
Method 2 is believed to be less accurate than Method 1, 
since it involves more approximations to the real gas 
properties. 

For other values of the Lewis Number, the effect of 
Lewis Number on the heat-transfer parameter was found 


by Method 2 to be best given by 


! 
l | 


Nu/+vV Re)/(Nu/vV Re), =1 = 


(L°*2 — 1) (hp “h,) (60) 
where the ‘‘dissociation enthalpy”’ /p is defined as 
hn = 3 c,.(—b/) = hat 2 e, (61) 


ATION 
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i.e., Ap is the dissociation enthalpy per unit mass of air 
in the external flow. The numerical results are plotted 
in Fig. 3 for comparision with Eq. (60). 

It was also possible to determine the effect of Lewis 
Number from the Method | computations by evaluat 
ing the additional term in Eq. (37) involving (L — 1) 
from the equilibrium properties of air. (The approxi 
mation for this evaluation is discussed in the Appendix.) 
Two such cases were computed, and the results are 
While 


there is some disagreement with the results of Method 2, 


plotted in Fig. 3 for comparison with Eq. (60). 


it is not too unreasonable considering the many different 
approximations involved in fitting curves to the func 
tions /, €,, ete. It is the authors’ opinion that the 
Method 2 solutions give a better indication of the Lewis 
Number effect for the equilibrium boundary layer, as 
embodied in Eq. (60). However, for Lewis Number 
unity the effect of pu variation is believed to be better 
given by the results of Method | [Eq. (5S) ], so that the 
total effect may be obtained by combination of Eqs. 


(58) and (60) in the form 


Nu V Re = 0.67( pus Pubtc)®* X 


11 + (L582 — 1)(hp hy} 2) 
The stagnation point heat-transfer rate for a 0.71 
thus becomes, by virtue of Eq. (45),” 
G = O.9A( picttic)®*"(psas)"* X 
11 + (2%? — 1)(hp hy){ (hs — hy) (du, dx) 63) 


It is interesting to note that the external flow proper 
ties are much more important than the wall values in 
determining the heat-transfer rate, so that the uncer 
tainty in the heat transfer is about 40 per cent of the un 


* For o not equal to 0.71, it is recommended that the factor 





0.94 be replaced by (0.76 07% 
1.0 +— 
08 + i | | | 
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Fic. 2. Correlation of the heat-transfer parameter Vu; V Re 
as a function of the pu ratio across the boundary layer, pu, Putty. 
for the equilibrium stagnation point boundary layer with L = 1. 
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certainty in the external viscosity. The physical reason 
for the importance of the external viscosity is that the 
growth of the boundary layer, and hence the heat trans- 
fer to the wall, depends mostly upon the external proper- 
ties. The analogy with turbulent boundary layers is 
easily seen. 

For a modified Newtonian flow, the stagnation point 
velocity gradient is 


(du, dx) (1 R)V2(p = Ps) / io (64) 


where FR is the nose radius and p,, is the ambient pres- 


sure. 


(9) THE ‘‘FROZEN’’ BOUNDARY LAYER 


When atomic gas phase recombination is negligible 
(t= 
reach the wall. 
recombination, the atom fraction at the wall will build 
up to the free-stream value. On the other hand, if the 
wall is cata/ytic, the atom concentration will be reduced 
to its equilibrium value at the wall temperature. Inter- 
mediate cases of wall catalyticity are of course possible, 


0), atoms diffusing from the free stream will 
If the wall is noncatalytic to surface 


but only these extremes were computed. 

Eqs. (34), (35), and (36) were solved with C; = 0 for 
various stagnation conditions and Lewis Numbers as 
discussed in the Appendix (Method 2). 

For Lewis Number unity, the effect of the pu varia- 
tion was very close to that found for the equilibrium 
boundary layer by Method 2, and could suitably be ex- 


pressed by Eq. (5S). For other values of L, the de- 
pendence could best be given by* 
(Nu /WRe)/(Nu/W Re), =; = 1 + 

(L968 — 1)(hp hy) (65) 


The calculated values are compared with Eq. (65) in 
Fig. 4. 

The difference in the exponent of Z for the frozen as 
compared with the equilibrium boundary layer [Eqs. 
(60) and (65)] is quite certain since exactly the same 
property variations were used in both cases, and also 
seems reasonable in view of the greater importance of 
diffusion throughout the whole of the frozen boundary 
layer. It can be seen, however, that for a Lewis Num- 
ber not too far from unity there is little difference in 
heat transfer for the frozen as opposed to the equi- 
librium boundary layer. 

A few cases for noncatalytic wall were also com- 
The resultant heat-transfer parameter could be 
the 


puted, 
given approximately by Eq. (62) with L = O—+e., 
heat transfer becomes proportional to h, — /p. 

A comparison of the distributions of enthalpy, tem- 
perature, and atom concentration for an equilibrium 
and a frozen boundary layer with identical catalytic 
wall and free-stream conditions is shown in Figs. 5 and 6. 
Both cases give very nearly the same heat transfer; 
however, the enthalpy distributions are slightly different 
(since 1 = 1.4) and the temperature and concentration 

* Note that Lees [see reference 7, Eq. (19)] suggested that, for 


the frozen boundary layer, the exponent of Z in Eq. (65) be 2/3. 
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distributions are markedly different, as is the mechanism 


of heat transfer. 


(10) FINITE RECOMBINATION RATI 

This most general case was solved using Method 2? 
with values of the recombination rate parameter (C 
varying from zero (frozen) to infinity (equilibrium), 
As for the frozen boundary layer, the wall may be either 
catalytic or noncatalytic, and both of these alternatives 
were calculated. It is to be expected, of course, that for 
large values of C; (near equilibrium) there should be 
little effect of wall catalysis, since few atoms reach the 
wall. 

The heat-transfer parameter for one flight condition 
and wall temperature is plotted in Fig. 7 for the com- 
plete range of C;. The solid lines are the total heat 
transfer for both catalytic (upper curve) and non- 
catalytic (lower curve) surfaces. For the catalytic wall, 
the fraction of heat transfer by conduction alone is 
shown by the dotted curve, so that the freezing of the 
boundary layer as recombination slows down (C; de- 
creasing) 1s easily evident. 

For either wall condition, C; must change by a factor 
of 10% in order for the boundary layer to change from 
substantially frozen to equilibrium throughout. Within 
this region of variation of C,, the boundary layer will be 
partly frozen (near the outer edge) and partly in equi- 
librium (near the wall). Since the recombination term 
[Eq. (56)] varies as 7~*, and the temperature changes 
by a factor of twenty between wall and external flow 
for the case considered, large variations in the recom- 
bination rate are possible across the boundary layer, 
thus permitting it to be partly frozen and partly in 
equilibrium. 

For the noncatalytic wall, the distributions of atom 
mass fraction for several recombination rate parameters 
are shown in Fig. 8S. For C, very large, no atoms reach 
the wall, all having recombined in the gas. 
values of Ci, some atoms reach the wall and, because 


For lower 


none recombine on the wall, a finite atom concentration 
builds up. For C, approaching zero, there is no re 
combination and hence no concentration gradients 
exist. 

It can be seen in Fig. 7 that a much lower value of C, 
is necessary to ‘“‘freeze’’ the boundary layer when a 
noncatalytic wall is used than would be the case other- 
wise. This is caused by the ‘damming up”’ of the atoms 
at the noncatalytic wall, resulting in greater recombina- 
tion because of high local concentrations. 

From Eqs. (57) and (63) it can be seen that, for a 
given flight velocity (hence 7°,), C,; varies as the square 
of the stagnation point density (and thus for strong 
shock waves, as the square of the ambient density), and 
also as the radius. Thus the boundary layer 


would become frozen at a high enough altitude, this 


nose 


altitude being less for small nose radii than for large. 
In order to change from a frozen to an equilibrium 
boundary layer, C; must change by 10‘, and thus the 
density by 10°, which is an altitude change of about 
100,000 ft. 
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(11) CONCLUSIONS 


The laminar stagnation point heat transfer in disso- 
ciated air can be given by Eqs. (63) and (65) for the 
equilibrium and frozen boundary layers, respectively. 
These results were computed for a Prandtl Number of 
0.71 and Lewis Number which 
throughout the boundary layer. 

The major deviation in the heat-transfer parameter 


for a was constant 


from the low temperature, perfect gas value is due to the 
The heat 
transfer [Eq. (63) ] is mainly dependent upon the value 


variation of pu across the boundary layer. 


of pu at the outer edge of the boundary laver. 

If the wall catalyzes atomic recombination, the total 
heat transfer is not much affected by a nonequilibrium 
state of the boundary layer if the Lewis Number is near 
unity. 

If the wall is noneatalytic, the heat transfer may be 
appreciably reduced when the boundary layer is frozen 
throughout—i.e., when the reaction 
time becomes much longer than the time for a particle 
to diffuse through the boundary layer. Since the ratio 
of these times depends upon altitude and nose radius, 
there is a scale effect which determines the chemical 


recombination 


state of the boundary laver. 


APPENDIX—DETAILS OF THE NUMERICAL SOLUTIONS 


Method 1 


The momentum equation for the stagnation point 
namely, 


boundary layer is given by Eq. (35) 


(1/2) ((p./p) — f,7] = 0 (A-1) 


fan)» 7 ug + 


If thermal diffusion is neglected (i.e., 1,’ = 0) then 
the energy equation in terms of the enthalpy becomes, 


from Eq. (37), 


[(/ o)gal, + Shy + 


1 a) Z[en(hi — hi) he\(Li — Vs, = 0 (A-2) 
or 
[(/o)(1 + d)g,], + fe, = 0 (A-3) 
where 
d = Dlc,,(h; — hi) hy |\(L;i — 1)(Os; Og), 
= (L — 1)D(h; — hj®)(Oc;/Oh), “( (A-4) 


where the subscript p denotes that the differentiation is 
at constant pressure and it is assumed 1; = L = con- 
stant for all species. 

In Method 1, which is suitable only for the equi- 
librium boundary layer, Eqs. (A-1) and (A-3) were 
solved simultaneously with the boundary conditions 


f(O) = 0, f,(0) = 0, 


f,,\ oO) = l 
g(O) = gy, g(o) = | 


The functions /, p, p and d were evaluated from the 
calculated equilibrium properties of air'! and by taking 
the viscosity to vary according to Sutherland's formula 


[see Section (5) above]. For given external (stagnation 
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point) flow conditions these quantities were plotted as 
functions of g. For numerical computation it was con 
venient to use analytic expressions of the following 


form: 
= Pl Pulte = (a V2) — (a2 g (A-5 
p;/p = 1 — y(1 — g) — ve] g (A-6 
d = (L — 1)3(h,; — h,®)(0c;/0h), = Bie” * (A-7 


The constants a, y, 6 in each expression were deter. 
mined by fitting these expressions to the equilibrium air 
calculations (see Fig. | 

Numerical solutions for this problem were obtained 
on an IBM 650 computer. The method of solution was 


to pick values of f,,(0) and g,(O) and integrate the 
equations directly, recording the resultant asymptotic 
values of f, and g for large values of 7. After three 
such integrations an interpolation will produce better 
values of /,,(0) and g,(0).. The interpolation procedure 
was repeated until the required conditions at “‘infinity 
were met—.e., f, > | and g— 1. This interpolation 
was made an integral part of the numerical program so 
by starting with three initial guesses for /,,(0) and g,(0 
the program would run automatically to completion. 


It should be noted that Eqs. (A-1) and (A-3) are 


formally identical with the stagnation point equations | 


solved by Cohen and Reshotko" except for the function 
d |Eq. (A-4)]. If ZL = 1, however, d = 0; thus by 
specifying 1 = 1,0 = 1,/ = 1 and p 

nation point solutions given by Cohen and Reshotko 
(A table giving the specific values 


p = g, the stag- 


could be duplicated. 


of the parameters for which solutions were obtained by | 


the method described above may be obtained directly 


from the authors. 


Method 2 

This method is a more general formulation in that it 
allows computation of the nonequilibrium boundary 
layer. As may be expected, however, it involves more 
approximations than the rather straightforward proce- 
dure of Method 1. 

In the nonequilibrium case, the concentration of the 
various species is not determined by the enthalpy and 
the (known) It is necessary, therefore, to 
add a continuity equation for each species. Further- 


pressure. 


more, it is convenient to express the thermodynamic 
properties in terms of the temperature and the con- 
centrations of the species. The energy equation should, 
therefore, be written in terms of the temperature. To 
make this problem tractable it was assumed that air is 
a diatomic gas composed of “‘air’’ molecules and ‘‘air’ 
atoms with properties properly averaged between oxy- 
gen and nitrogen. The dissociation energy of an air 
atom was taken to be the average dissociation energy in 
the external flow [see Eq. 59)]. With this assumption 
the problem is reduced to the simultaneous solution of 
three equations (momentum, energy and atom con- 
centration) and the thermodynamic properties are to be 
expressed in terms of the temperature and atom concen- 


tration. 
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THEORY OF STAGNATION 
The momentum equation is still Eq. (A-3), with 
p./p = (1 + ca)/(1 + cas) (A-8) 
Pu Pickte = [1/11 + ca) ]° (0, 0) FO) (A-9) 
wher¢ 
Fi6 (7,0/300)*'* [413 (7,8 + 113)] + 


3.7(7.6/10,000)? — 2.35(7.6,/ 10,000) 4 


and the stagnation temperature 7), is given in degrees 
Kelvin. The function F(@) is a fitted curve giving the 
temperature dependence of the viscosity under the as- 
sumption that the atoms and molecules have the same 
collision cross-sections. 

The energy equation in terms of the temperature is, 
irom Eq. (36) with L,’ = 0, 


cl'a)0,|, + cf0, + (1/c)6, X 
V(Cpi Ene) Le CisSin + 4 2(du,/ dx) : x 
V(w;/p)((h® — hy)/EpeTs.| = 0 (A-10) 
where ri Ca/ Cos 


With the assumption of a simple diatomic gas and taking 
L; = L = constant, the third term may be rewritten as 


CpM), Cpu 


(L/ 7)0,C4y|(Cpa 

and the fourth term, using Eqs. (56), (57), and (61), 
and taking h4 = hy,, becomes 

Ci(hp ‘ee & ) [(c4* — Cap") /6°5 

Cow © (R/M)} (7/2) + ce "| 


Ni Ww 


where the exponential is the vibrational heat capacity 
and 7, ~ SOO°K. for air; also 


Cpa = (0 2)[R/(.W/2)] 
Hence, 
C= 6, Eye = (10/7)ca + 1 + (2/77 NL = 4) 
(A-11) 
CypA — Cym) Cou = (3/7) — (2/)e es (A-12) 


For computation then the energy equation becomes 


[(cl/o)0,], + cf0, + e(Ll/c)6,c4, + CiCom = 0 (A-13) 
where 
l pit Pht (1/(1 + c4)]*/7(6,./0) F(0), 
see Eq. (A-9) 

é (Cyd — Com) /Cpw = (3/7) — (2 ieee 
C == C/€yn = (10/ tea + 

51 + (2/7)e Or F(A — ca) 
C, parameter, see Eq. (57) 
C; hp Enel, = parameter, see Eq. (61 
m (Ca> — Car”) /0**(1 + c4) 


The equilibrium atom mass fraction c,4z can be deter- 
mined from reference 11. For computation c4, was 
approximated by 


(A-14) 


where C3 is a constant. 


POINT HEAT TRANSFER 

The continuity equation for atoms was written im 
terms of the atom mass fraction C4 instead of the nor- 
Thus Eq. (34 


malized atom mass fraction s. becomes 


or 


lL ove =o cat. lm Cm = U A-16) 
Method 2 for the nonequilibrium boundary layer is 
the simultaneous solution of Eqs. (A-1), (A-13) and 


(A-16) with the boundary conditions 
AO) «= & 0, f(a ] 
60) = Oy, 


f(Q) = 
@( ao) = |] 


c4(O) = O for catalytic wall ( 
C4,(0) = 0 for noncatalytic wall} 


Solutions were obtained on a digital computer using an 
iterative procedure similar to that used in Method 1. 

The limiting case of the equilibrium boundary layer 
was obtained by Method 2 by eliminating the term Cym 
between Eqs. (A-13) and (A-16) and solving the result 
ing equation simultaneously with Eq. (A-1), taking the 
equilibrium atom concentration as a known quantity in 
the form of Eq. (A-14). The limiting case of the frozen 


boundary was obtained by putting C; = 0.” 
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SUMMARY 


Phe results of an experimental investigation of the laminar heat 
transfer at the stagnation point of a blunt body in partially dis 
sociated air are presented and are compared to the theoretical 
Riddell 


sented for air temperatures as high as 8,000°K 


treatment of Fay and Heat-transfer results are pre 
where more than 
50 per cent of the molecules are in the dissociated state. The ex- 
periments were performed in a shock tube and the new experi 
mental techniques and principles are discussed briefly Simula 
tion of flight stagnation conditions at velocities up to satellite ve 
locity of 26,000 ft 


data has been obtained over a large altitude range at these veloci 


sec. is shown to be possible in shock tubes and 


ties 
SYMBOLS 
A = ratio of thermal conductivity and diffusivity 
K i/Km) WV ku /k 
a = speed of sound 
C, = heat capacity 
D = diffusion coeflicient 
E = electrical potential 
h = enthalpy per unit mass 
hy = average dissociation energy per unit mass of atoms times 
the atom mass fraction in the external flow 
I = electric current 
AK = thermal conductivity 
k = diffusivity K/pC 
I = Lewis Number DpC,/K 
l linear dimension, thickness 
VW = Mach Number 
Vu = Nusselt Number qx/K AT 
p = pressure 
q = specific heat flux 
R = resistance 
Re = Reynolds Number pul) 
1 = temperature 
t = time 
u = velocity 
x = linear dimension 
a = temperature coefficient of resistivity 
> = ratio of specific heats 
u = viscosity 
p = density 
a = Prandtl Number Cyp/K 
T = characteristic time of experiment 
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Subscripts 


27 = insulator 
m metal 
o = constant or steady state 


s = stagnation conditions 
= wall 
1 = initial conditions 


2) = conditions behind moving shock 


= driver conditions 


INTRODUCTION 


TT" PROBLEM of aerodynamic heat transfer at the 
stagnation point of a blunt body in hypersonic 
flow, characterized by large temperature gradients in 
the boundary layer has received a considerable amount 
of attention in the literature (see, for example, Lees 
Cohen Reshotko,*? Sibulkin,' 
cently, in reference 3, Fay and 


) 


and others). Re 
Riddell 
including the 


and 
have con 
sidered the more general problem by 
effects of variable fluid properties in the boundary 
layer, diffusion of atoms, and chemical processes such 
as dissociation and recombination. In this treatment 
the heat transfer at the stagnation point has been solved 
for several cases; (1) a boundary layer in thermochemi- 
cal equilibrium throughout, (2) a boundary layer whos¢ 
composition is “‘frozen’’ at the atom and ion concen- 
trations that exist at the edge of the boundary layer 
and (3) a boundary layer wherein the concentration of 
each chemical species is permitted to adjust as dictated 
This 


presents the methods and results of the experimental 


by a finite atomic recombination rate. paper 
attack on this problem. 

At the present time, much effort is being spent on the 
development of a facility which is capable of repro- 
ducing, in the laboratory, the environmental conditions 
encountered by a vehicle traveling at hypersonic ve- 
locities. Some of the schemes generally under consider- 
ation are hypervelocity guns, wind-tunnel modifica- 
tions, free flight missiles and several versions of the 
shock tube. Of these, only the shock tube has the 
ability to produce the correct conditions of temper- 
ature and pressure while holding the model stationary 
This 


method of operation has the advantage of simplifying 


i.e., by moving the air as is done in wind tunnels. 


the instrumentation, which is very difficult in any hyper- 


t The “frozen” and “equilibrium” situation can best be dis- 


tinguished by considering the relaxation time of the chemical 
as compared 


process—i.e., recombination of dissociated atoms 


to a characteristic flow time, such as the time required for a 


particle to flow through the boundary layer. The frozen flow is 


defined as TfJowK Tree. 


mb, and equilibrium as tyjow >> Tp 
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ress of the wave and the 

The gases which were originally separated 

separated by the contact discontinuity 

distribution at several typical times 
and ¢; from (b 


expansion 


by the 
interlace < 
in the shock tube 


sonic facility. In addition, the low cost of operating a 
shock tube makes it a very attractive instrument for con- 
ducting basic high temperature gas dynamic research. 
The shock tube heat-transfer experiments described in 
this paper involved the development of new instru- 
mentation capable of microsecond response. For sur- 
face temperature measurements, the thin film resistance 
thermometer, as used by other investigators for similar 
purposes, but at low Mach Numbers, was developed 
for use under the extreme temperatures and pressures 
encountered at the stagnation point at very high Mach 


Numbers. * 


In addition, a new type of meter has been 
Heat- 


transfer measurements have been performed over a 


developed based on the calorimeter principle. 


wide range of shock tube pressures and temperatures, 
simulating hypersonic flight at velocities up to about 
26,000 
theoretical analysis of Fay and Riddell*® and indicate 
with the 


ft./sec. The results are compared with the 


reasonably good agreement equilibrium 


boundary-layer solutions. 


THE SHOCK TUBE 


The shock tube, as used for these experiments, has 


been described in the literature.‘*~* Schematically, the 


shock tube is shown in Fig. 1. A strong shock wave is 


created by the sudden removal (rupture) of a dia- 


phragm separating a high- and low-pressure gas. 


Compression waves advance into the low-pressure gas 
and quickly (in a few diameters) steepen to form a 


shock wave. Simultaneously, the high-pressure region 


\s well as electrical-conductivities of air that are comparable 


1 


to the resistivity of the gage resistance element 


also expands, but the interface between these two gases 
moves at the particle velocity which is slightly slower 
than the shock 
shock wave runs away from this interface, or contact 


wave velocity. Consequently, the 
discontinuity, creating a continually growing region 


When the 


shock wave advances into the low-pressure gas (air for 


as the shock wave moves down the tube. 


the present experiments) it behaves according to the 
laws of conservation of momentum, energy, and mass, 
resulting in a gas which has been compressed, heated, 
and set into motion at almost the velocity of the shock 
wave. By a change of coordinate system —1.e., by con 
sidering the problem from the point of view of an ob 
server traveling at the shock wave velocity -this 
process is reduced to a more generally. known form ot 
the Rankine-Hugoniot equations. 

The conventional method of presenting the processes 
occurring in a shock tube is the x-/, space-time, diagram 
in the center part of Fig. 1. The diaphragm breaks 
at distance and time equal to zero. In practice, the 
interface is found to spread, probably due to turbulent 
mixing, and is indicated as a region rather than a link 
interval be 


Downstream in the shock tube, the time 


tween the passage of the shock wave and the arrival 
of the mixing region gives the experimenter a region 
which has been uniformly heated and compressed in 
which to make measurements. 

The pressure at every point in the tube at several 


typical times is presented in the lower part of Fig. 1 

It is possible to produce a wide range of conditions 
of temperature, pressure, and velocity in shock tubes 
limited only by the strength of the shock wave that 
can be created in a given gas. The strength of this 
shock wave can be expressed as a function of only the 
pressure ratio across the diaphragm and the ratio of 
the speed of sound in the driver and in the driven gas, 
as well as the specific heat ratio in the driver gas, 4, as 


follows: 
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Fig. 2 shows the theoretical strength of the shock 
waves produced in air driven by cold hydrogen and also 
by a mixture of oxygen, hydrogen and helium in which 
an increase in the speed of sound ratio is achieved by 
combustion. The experimentally achieved shock wave 
strengths agree quite well. 

The equipment used for these experiments is capa- 
ble of producing shock waves traveling at over twenty 
times the speed of sound in room temperature air. 
If the air were an ideal diatomic gas (y = 1.4) this 
strength shock wave would heat the air to a temper- 
ature of about 24,000°K. 


goes into dissociating the oxygen and nitrogen mole- 


However, energy actually 
cules rather than into raising the temperature. Using 
the properties of air in thermochemical equilibrium, 
either from the National Bureau of Standards tables,’ 
or from the data calculated by the Rand Corporation,’ 
in the solution for the conditions behind the shock 
wave, the air is actually found to attain a temperature 
of only 9,000°K. 

Correlation of shock 
tubes with quantities calculated from theories requires 
a precise knowledge of the state of the gas in the shock 


physical measurements in 


tube at the point where the measurements are being 
conducted. For this purpose the initial pressure in the 
driven gas and speed of the moving shock wave are 
measured. From these two measurements, the proper- 
ties of the gas at any part of the shock tube at any time 
can be calculated. For this purpose, the Avco Research 
Laboratory has constructed a Mollier diagram for air® 
from the National Bureau of Standards tables. Refer- 
ence 9 contains curves of the solutions of the Rankine- 
Hugoniot equations for equilibrium air for typical 
shock tube and hypersonic flight conditions. 

It must be pointed out that in the example discussed 
above, the speed of sound in the gas heated by the in- 
cident shock wave is approximately 6000 ft. per sec. 
due to its high static temperature, so that the air even 
though traveling at almost 20,000 ft. per sec. has a 
flow Mach Number of only about three. 


SIMULATION 


At the stagnation point of a blunt body in hypersonic 
fight, the usual separation of theoretical aerodynamics 
into a viscous and inviscid problem is exactly applicable 
to the experimental situation. As long as there are no 
interactions between the shock wave and the boundary 
layer, such as when the boundary layer is of the same 
thickness as the detachment distance, 
wave boundary-layer intersections are present or when 
the boundary-layer 


where shock- 


flow separations are involved, 
problem can be considered experimentally apart from 
the inviscid flow field and Reynolds Number simulation 
does not have to be achieved. 

The temperature, pressure, and velocity environ- 
ment created in the shock tube can be applied directly 
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to the solution of the stagnation point heat-transfer 


problem. Aerodynamic heat transfer at the stagna- 


tion point, by conduction, convection, and diffusion of | 


air molecules and atoms is simulated in a scale test if 
the conditions at the edge of the boundary layer in the 
vicinity of the stagnation point, including the pressure 
distribution and the wall conditions, are correctly re- 
produced. 

It can be shown that the heat-transfer rate through 
a laminar boundary layer at the stagnation point at 
low speeds, where the product pu in the boundary layer 


does not vary from the free stream to the wall and the 
. | 


density varies inversely with the temperature,” can be | 


expressed in a nondimensional form as 
Nu ~~ Re = constant = 0.763093 (2) 


The numerical solution of Fay and Riddell* for the 
equilibrium boundary layer with variable fluid proper- 
ties can be correlated by considering the pu ratio across 


the boundary laver, by the expression 


Nu/V Re = 0.7630% 4; [(pu)s/(pu)w]°4 X 
[1 + (195? — 1) (hp/hs)]j (3) 


The heat transfer at the stagnation point can be 
written from this expression, assuming the pressure 
distribution can be expressed by the Newtonian as- 


sumption, as 


L) hn / Rd 
VR) (2p. ps)"/44 (4) 


gq = 0.763) (pu) u?*'(pu)es4 [1 + (195? - 
[(hy — ly) /o%®] 


This equation makes the heat-transfer simulation at 
the stagnation point dependent on the inviscid flow 
parameters, the wall conditions, and the pressure gra- 
Note that the wall conditions enter as (py),,!:! 
— h,) and affect Eq. (4) very weakly, so that 


dient. 
and (h 
high cooling, as long as the cooling is large, is sufficient. 
In addition, the transport properties of air —i.e., u, A, 
and D —enter the solution in the form of the 
Number, L = DpC,/K, Number, 
Cpu /K. 

The inviscid flow conditions produced at the stag- 


Lewis 


CC = 


and Prandtl 


nation point of a model in the shock tube are used to 
simulate the hypersonic flight situation. However, 
the stagnation enthalpy is achieved in a manner dif- 
ferent from flight -1.e., in a two-step process: first, 
through the conversion of the kinetic energy of the 
incident shock velocity, “,° 2, into the total enthalpy of 
the hot gas region (as considered in the shock fixed 
coordinate system) and, second, through the conver- 
sion of the total enthalpy of the moving hot gas 
into the stagnation enthalpy in the shock tube fixed 
coordinate system. The initial pressure level of the ex- 
periment can be adjusted to obtain the proper stagna- 
tion density, p;, in the inviscid flow. 

The boundary layer on a model in the shock tube 
under the same stagnation conditions as in hypersonic 
flight will be an exact geometric scale reproduction of 
the flight case if the two situations are chemically 


similar. This means that the flight and shock tube 
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Schlieren photograph of the flow over the nose of a 


FIG. 5 
kemisphere-cylinder model in a 1'/2 in. inside diameter shock 
tube Note fine Mach lines whose inclination to the body is 
accurately measurable 
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boundary layers must be either both in thermochemical 
equilibrium throughout or both completely “frozen.” 
Between these two extremes the scaling of the recom 
bination parameter discussed by Fay and Riddell* is 
applicable. The implication of the latter type of simu- 
lation is that the stagnation enthalpy and the atom 
diffusion process through the boundary layer are cor- 
rectly duplicated in the experiment. The only other 
variable which must be determined to complete the 
simulation principle is the velocity gradient at the 
stagnation point or pressure distribution. 

For blunt bodies in hypersonic flow, the Newtonian 
approximation predicts that the pressure distribution is 
independent of Mach Number. This approximation 
to the pressure distribution has been shown experi- 
mentally to be correct even at surprisingly low Mach 
Numbers for very blunt bodies. It remains only to be 
shown that this approximation can be made in the 
shock tube. Because of the difficulty of making reli- 
able pressure measurements in the time scale of the 
shock tube an optical scheme from which the pressure 
distribution can be calculated has been developed 

Fig. 3 shows a schlieren photograph of the flow over 
a hemisphere-cylinder model in a 1! 2 in. inside diameter 
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circular shock tube. The surface of the model was 
roughened slightly so that weak disturbances—-.e., 
Mach waves—-propagate outward from the body. 
By following an isentrope on a Mollier diagram from 
the stagnation point, the pressure at every point 


where a Mach angle is measured can be calculated for 
equilibrium air.* 

Fig. 4 shows the pressure distribution deduced in 
this manner from an experiment performed at shock 
Mach Number 8.8 in an initial pressure of 10 cm. of 
Hg of air (flow Mach Number about 2.3). The pressure 
distribution does not the Newtonian estimate 


exactly. However, the sonic point, as indicated from 


follow 


the first visible Mach line, seems to occur very close to 
the predicted point. In addition, the slope of the pres- 
sure distribution at this point is very close to the esti- 
mated slope, leaving little doubt that the pressure dis 
tribution inside the sonic point is close to Newtonian. 

The fact that the pressure distribution in the vicinity 
of the stagnation point on a hemisphere cylinder model 
can be approximated by the Newtonian estimate, even 
for relatively weak shocks, makes the shock tube a 
very useful tool for experimental boundary-layer stud- 
ies. Simulation of stagnation point heat transfer in a 
shock tube thus becomes completely independent of 
flow Mach Number and simply a matter of geometric 
sealing and duplicating enthalpy and pressure. 

Fig. 5 shows the stagnation conditions attainable 


in a shock tube as compared to those encountered in 


rhis calculation depends on how the speed of sound is de- 


, equilibrium or ‘‘frozen”’ 


3 and 4, 


fined along the isentropic expansion—41.« 
it stagnation point dissociation. For the example of Figs 
the degree of dissociation is negligible so that variances with these 
definitions are eliminated. The fluid properties in an isentropic 
expansion in equilibrium air can be calculated from the use of an 
19 for the case 
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FLIGHT VELOCITY -FT/SEC x 1073 


Shock tube stagnation conditions compared to the 


Fic. 5 
stagnation conditions encountered in hypersonic flight. In ad 
dition, conditions attainable by a 36:1 area ratio nozzle are 
indicated for one initial pressure. Typical trajectories'® of 


several types of missiles are also shown 
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Fic. 6 Nomograph * for calculating the time interval between 
the arrival of the shock wave and the driver gas ata given 
station in a shock tube. Calculations have been carried out 
considering air an ideal gas = 1.4), and in equilibrium For 


equilibrium air the testing time is greatly reduced 
hypersonic flight, on a pressure-temperature grid. The 
conditions attainable in a straight shock tube without 
the addition of expansion nozzles are plotted. Trajec- 
tories of various types of proposed hypersonic vehicles, 
as given by Masson and Gazeley,'® are shown to fall 
well into the accessible range. In particular, trajec- 
tories of ballistic bodies fall into the range covered by 
the straight shock tube. For better simulation of the 
flight of glide vehicles, satellites and continuous flight, 
the shock tube gases can be expanded in a nozzle ex- 
tending the shock tube capabilities to somewhat higher 
altitudes. !! 

Most of the experiments reported in this paper were 
performed in 1! , in. inside diameter shock tubes, vary 
ing from 15 to 30 ft. in length. 


MEASUREMENT TECHNIQUES 


Because of the extremely short steady-state times 
available in shock tube experiments, instrumentation 
presents a problem. Fig. 6 gives the testing time in 


shock tubes under typical conditions. The departure 


* Enter abscissa at shock Mach Number and find intersection 
with initial pressure curve. Locate values of left-hand ordinates, 
it the latter intersection Line up the value of 7/L with the 
length of the shock tube on the right-hand ordinate and read the 
testing time 7 on the center ordinate. Experimental values of 
7 are found to be approximately 40 to 70 per cent of the values 


calculated this way 
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from the ideal gas gives shorter times due to the higher 
real gas density ratios. In practice, it is found that 
only 40 to 70 per cent of the theoretical time is actually 
available, probably due to mixing at the interface and 
to attenuation of the shock strength. 

Typical of the instrumentation problem is the meas- 
urement of the velocity of the incident shock wave, 
from which the stagnation conditions are calculated. | 
In this experiment, the passage of the shock wave at a 
number of points along the length of the tube is re. 
corded by ionization-sensitive pickups. The signals 
are displayed on the face of an oscilloscope on a sweep | 
that has been extended to cover very long times, yet 
is readable to very short intervals. An example of the 
trace of this measurement is shown in Fig. 7. Accu- 
racy better than +2 per cent for each interval can be 
obtained, depending on the accuracy of the signal and 
pickup spacing. 

In order to make significant heat-transfer measure 
ments in shock tubes an instrument capable of micro- 
second response had to be developed. The most prom 
ising schemes considered were of two types, (a surface 
temperature history measurements, and (b)  calori- 
metric measurements. 

Secause of the difficulty of shock tube instrumenta 
tion, it was felt highly desirable to make several inde- | 
pendent measurements of the heat-transfer phenomena | 
Variations of the resistance thermometer principle were 
found satisfactory for both above schemes. 

If one considers the one-dimensional heat conduction 
problem for a finite thickness of a metal, backed up 
by an infinitely thick insulator, with a constant heat 
flux, g,, entering the metal surface, then the temper- 
ature response of the interface between the metal and 


the insulator can be written as follows: 














. o 
T Al. t) Ae. £) am 
2V Rk, / = 1 “ 
K (1+A)m=™ \1+A y 
as i 2n + | . (2" + 1) 7 
T V k, 2V &. 


If the metal thickness, /, is very small, that is, very 
much smaller than the characteristic diffusion length 
in the insulator, / V k,7,7 then the temperature re 
sponse of the metal essentially is the same as the re 


+ And therefore also much smaller than the same length for the 


metal, / V k,»7, which is larger than /;. 
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Ss, yet 
of the 
Accu- 
‘an be 
il and Fic. 7. Typical oscillograph from the shock velocity timing 
ircuit. The oscilloscope sweep starts at the top and makes a 
mplete cycle, left to right and back, every 100 microsec. The 
isure- bright spots are produced by the output from a crystal oscillator, 
Sa spaced exactly 10 microsec. apart. The time of arrival of the 
nicro- : ae 
shock wave at each successive ionization-sensitive pickup is 
prom- lenoted by the sudden break in the sweep (change in gap resist 
orfore nce caused by the ionized gas Sample oscillograph has signals 
ait rom 12 pickups. The pickups are spaced approximately a foot 
‘alori- part 
lenta > . £ * 
sponse of the interface temperature,* and can be re- 
inde- , 
| duced to 
nena 
were 1, (1, t (2V ki V rK3)Q.V ft (6) 
ti For the case of the thin resistance thermometer, a 
ction , ‘ . 
; current is passed through a thin metal film through 
d up ' ; 2 : = 
; ; which the heat is transferred to the insulator. The 
1€a P . ° 
resistance change of the metal, interpreted as the sur 
1per- . , ; 
; | face temperature of the insulator, yields the heat 
am . : . ; : - 
transfer rate from Eq. (6). In terms of measured 
quantities, the heat-transfer rate can be written as 
do = [1 /k”UIR)] (E/Vt) (7) 
Kim ‘ / ? 
where k” = 2a,Vki/Vak (S) 


represents the physical properties of the resistance 
thermometer. 
The other extreme compared to the very thin film 
situation is a metal layer whose thickness is larger than 
| the characteristic diffusion length of the metal—i.e., 
> Vkn»r. Such a gage is one in which essentially 
all the heat stays in the metal—.e., a calorimeter. By 
expanding the expression for the temperature and re- 
sistivity of a thin metal lamina and integrating over 
the metal thickness, an average temperature can be 
found from which the response of the resistance element 


= to the heat flux is described as follows: 


(,)) 
Jo = [1/(R)] (pCl/ a), (dE dt) (9) 
Cry ities - , 
tl Chis expression is good as long as there are no two- 
gth ; ? ; ‘ ; ‘ . 
2 dimensional effects——that is, each lamina is at a uni- 
re ¢ 
form temperature. 
re 
rhe error introduced by the heat capacity of the thin metal 
the film has been considered in reference 19 and is negligible for the 


geometry used in these experiments. 


Both of the principles described above can be ex 
tended to the case of nonuniform heat flux in time—i.e., 
g = g(t). The calorimeter integrates the total heat 
input into the metal. Therefore, the heat-transfer 
rate at any time is given by the instantaneous slope of 
the voltage signal with time. 

For the thin resistance thermometer, the surface 
temperature response to a heat input arbitrarily chang 


ing with time can be written as 


Vv k, G\T) = 
i (oe, 2) = Tt) = dr (10) 
KV r- V 


This expression can readily be inverted to solve for 
the time-dependent heat-transfer rate as a function of 
surface temperature, which, in terms of the temperature 
coefficient of resistivity and the measured quantities, 


vields 


2K * E(t) — E(r) 
a(t) j t) li | _ 


R"(IR)e i Vt. 


This equation has been set up for solution by digital 
computer for reduction of the data obtained from the 
photographic oscillograph records. 

The thin resistance thermometer as used by Black 
man'? and others!®'4was developed mainly as a trigger 
device. At the stagnation point of models under the 
extreme conditions created by strong moving shock 
waves considerable development was required before 
quantitative data was obtained. It can be shown that 
a thin gage can be up to one micron thick, considering 
times of the order of 100 microsec. The metal films 
produced were measured by interference microscope to 
be approximately 0.3 microns thick. 

The production of these thin films required much 
development. Evaporated films were found to be 
unsatisfactory because they did not adhere to the 
insulator sufficiently. The most promising results 
were obtained through the use of a solution of organic 
platinum and gold compounds of resinous character in 
oils and other solvents. This solution, manufactured 
by the Hanovia Chemical and Manufacturing Com 
pany, is called Liquid Bright Platinum No. 05. The 
solution is brush applied in the desired geometry and 
fired in a ventilated oven to about 1,200°F. for a few 
minutes. In order to produce highly erosion-resistant 
films, several coats, each individually fired, were re 
quired. 

In the analysis of the data from thin heat-transfer 
gages, the characteristics of the gage are represented by 
the diffusivity and the thermal conductivity of the in 
sulator and the temperature coefficient of resistivity of 
the metal. The latter, as measured in a laboratory 
oven from 0° up to 200°C., when combined with the 
manufacturer's data for the properties of the insulator 
i.e., pyrex——gives a value of k”, defined by Eq. (8), of 
0.00765. Because the metal fuses into the outer layers 
of the insulator, the thermal conductivity of the pyrex 
could be greatly in doubt. To eliminate the uncer- 
tainty in the value of the calibration constant, k”, due 
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Fic. 8. Sample heat-transfer gages used. The thin gage is 
on the right and the calorimeter gage on the left. Both gages 
are mounted on ground and polished pyrex hemisphere-cylinder 
models. Gages normally plug into sting model support 
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Fic. 9. Heat loss from calorimeter to the insulator for various 
heat-transfer gage materials. Calculations are made assuming 
a perfect contact between the metal and insulator. 
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to this effect, a method of dynamically calibrating the 
gage in the true scale of the shock tube experiment was 
developed. 

When a pulse of electrical energy is passed through 
a resistance thermometer, the resistance of the element 
will change in proportion to the power dissipated. If 
this pulse is made comparable in energy and time dur- 
ation to thermal energy transferred during a shock 
tube experiment, the response of the resistance ther- 
mometer is directly comparable. As a result, from the 
known energy in the pulse, the properties represented 
by the value of k” can be calculated. The average 
value of the calibration constant, determined from the 
electrical calibration, brackets the calculated value 
within plus or minus 15 per cent. 

It was found to be a very difficult task to develop thin 
resistance thermometers which could withstand the 
extreme conditions, such as erosion due to dust and 
diaphragm fragments, without serious damage, as evi- 
denced by a change of resistance. As a result, gages 
are calibrated and used for one run mostly, except 
in the low shock strength and low density region. There 
is considerable evidence that the damage which does 
occur happens after the interface arrives, that is, after 
the data has been taken. Typical heat-transfer gages, 
both the thin and the calorimeter type, are shown in 
Fig. 8. 

Difficulties encountered with the thin resistance ther- 
mometers,* particularly at high Mach Numbers, led to 
the development of the calorimeter gage. Considering 
Eqs. (5) and (9), it can be determined that a metal 
of one characteristic thickness will absorb a large 
fraction of the heat flux incident upon it in 100 microsec. 
For platinum this thickness is approximately 0.003 cm. 
Such a strip of metal gives a response of one millivolt to 
a heat input of one kilowatt per square centimeter, 
which is the order of the minimum heat-transfer rates 
of interest in these experiments. 

Because the gage is not infinitely thick, there will be 
some heat loss to the back-up insulator. This loss can 
be calculated and corrected, assuming a perfect contact 
between the metal and insulator. This loss, as a fune- 


tion of gage thickness and testing time is shown in 


Fig. 9. In an actual gage there is a small air space 
between the two materials, which materially reduce this 
loss. 

The manufacture of calorimeter gages has been 


simplified to make the gages expendable, an almost 
necessary condition when making experiments with 
strong shock waves. A die-cut piece of platinum is 
cupped and mounted on a ground hemisphere-cylinder 
pyrex model by soldering it onto areas previously 
painted with the same platinum used for the thin gages. 
A four lead pattern, two to introduce the current and 
two to measure the resistance of a specific area, has 
been used to minimize the edge effects and aerodynamic 
disturbances caused by the wires and soldering. An 


* Such as the large amount of work involved in the reduction 


of the oscillographic trace to heat-transfer data 
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does Fic. 10 Calorimeter heat-transfer gage i exploded view 
after ind (b) schematic. The gage area, between the two thin poten- 

f tial leads, covers an area of +20° surrounding the stagnation 
Pe re . . Ps 
ARES, point. Heat-transfer data is corrected for the average over the 
vn in | gage area based on an assumed distribution (Newtonian pres- 
ther- | 
ed to , : 

i exploded view of a typical heat-transfer gage and cir- 
ering | ° 5 i aie ha 
i cuit schematic are shown in Fig. 10. The gage covers 
neta = . : : . - 
| | a finite surface surrounding the stagnation point. The 
arge ° ° ° 
- average heat-transfer rate over this surface has been 
psec. - . ‘ 
calculated to be 94 per cent of the stagnation-point 
em. ter 
It rate, assuming a Newtonian pressure distribution. A 
it to ° : ‘ 
corresponding correction 1s made. 
eter, ; , : 
ates lhe properties of the gage material and insulator of 
interest for the calculation of heat-transfer data from 
. “1 ter o -re > are C1ven ji ) ( 
il be calorimeter gage response are given in Eq. (9). 
can | For the thicknesses of metal used, bulk properties 
tact | are considered reliable. For chemically pure plati- 
inc- num, it has been estimated that the important proper- 
1 in ties, i.e., the heat capacity, density and temperature 
yace coeflicient of resistivity are known within an accuracy of 
this plus or minus 2 per cent. 

The electrical calibration, previously described for 
een | the thin heat-transfer gages, has been successiully used 
lost ; to check the properties used for calorimeter gage data. 
7 , ee 
ith | rhe low resistance of the thermometer element, ap- 
11S | proximately + X 10~* ohms, dictates the use of a 1,000- 
der | amp. pulse to produce energy dissipation in the gage 
sly comparable to the shock tube test transfer. It is 
yes. difficult to produce a reasonably square pulse of this 
ind current amplitude, but a pulse whose current variations 
1as are known can be used to obtain the material properties. 
mc The properties thus determined check the manufac- 
An turer’s data within plus or minus 10 per cent. 

The details of analysis and experiences collected in 
ion the development of shock tube heat-transfer techniques 


are reported in reference 21. 
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RESULTS 


rhe stagnation point of a blunt body has advantages 
which make it attractive for heat-transfer experiments. 
First, the most severe laminar heating occurs at this 
point; second, direct simulation of flight conditions 
i.e., correct stagnation pressure and enthalpy inde 
pendent of the Mach Number of the flow surrounding 
the model—can be accomplished at this point; and 
third, the problem of conductive, convective, and dif 
fusive* heat transfer can be treated theoretically, 
limited only by the uncertainty of the transport 
properties at very high temperatures. However, cal 
culations have been carried out over a range of vis 
cosity, Prandtl Numbers and Lewis Numbers large 
enough to cover the most reasonable estimates that 
have been made for these parameters. This uncer- 
tainty in the calculated heat transfer is of the same order 
as the experimental scatter. 

In addition, at the stagnation point of small models 
the radiative heat transfer is minimized by the small 


scale, whereas the aerodynamic heat transfer increases 


* Heat transfer due to diffusion of atoms 
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Fic. 11 Typical oscillographs ind data On all three oscil 


lographs time increases from left to right ‘he sweep speed is 


10 usec. /interval The oscilloscope is triggered before the mov- 
ing shock arrives at the stagnation point The break in the 
horizontal trace denotes the arrival of the shock at the gage 
The break of the curve accompanied by the end of the constant 
heat transfer denotes the arrival of the interface a) Record 
from a thin heat-transfer gage at a relatively low shock Mach 
Number b) Record from a thin heat-transfer gage to a higher 
shock Mach Number. Note that the extreme conditions at the 
high stagnation enthalpy degrade the quality of the signal c 


Calorimeter gage record 
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inversely with the radius ratio of the scaling. Thus, for 
a small, say 1/2-in., diameter hemisphere model com- 
pared to reasonable nose sizes, the ratio of radiative to 
aerodynamic heat transfer is suppressed by a very large 
factor, leaving measurements of purely aerodynamic 
heat transfer at high stagnation temperatures, even at 
11,000°K. 

The experiments were performed at the stagnation 
point of 1/2 
with either thin 
mounted at the stagnation point. 
were performed on 1-in. diameter models in a 
100 ft. long shock The 
ground and optically polished to the proper nose radius. 

Oscillograms of the heat-transfer 
gages during a typical experiment are shown in Fig. 11. 


in. diameter hemisphere-cylinder models, 


or calorimeter heat-transfer gages 
A few experiments 
f-in. 
models were 


tube. pyrex 


response from 
The thin gage has a parabolic response with time after 
the shock wave arrives at the gage. The oscilloscope 
is triggered a few microseconds before the arrival of the 
shock wave at the gage in order to record the initial 
portion of the curve. The output of the resistance 
element can be plotted against the square root of time, 


yielding directly, by use of the result of the gage cali- 


SEGAL 


L9YSS 


EBRI 


ARY, 


SCIENCES Fk 


bration, the unknown heat-transfer rate, 4g. 
records can also be analyzed by reading the output at 
many time intervals, and solving Eq. (11) for the heat 
transfer rate, g(/), as a function of time. The result of 
such an analysis is plotted in the left half of Fig. 11 
An oscillogram from a typical calorimeter gage js 
shown in the lower half of Fig. 11. 


variation of the output with time yields (from the values 


Che | 


The instantaneous | 


of the physical properties of the platinum strip as | 


supplied by the manufacturer) the unknown heat 
transfer rate, g(t), at any time. 

Fig. 12 presents the results of the experiments that 
have been performed. Three series of experiments at 
convenient initial shock tube pressures simulate hyper 
sonic flight at altitudes of roughly 20,000, 70,000, and 
120,000 ft. 
below the onset of chemical effects, at stagnation tem 


The velocity range covered extends from 


peratures of about 3,000°K., 
required by a satellite to counteract the earth’s gravi- 
field. 


calorimeter gages. 


tational The data were mostly taken with 


thin heat-transfer gage experiments. There seems to be | 


a reasonable agreement between data from the two 
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Measured stagnation point heat-transfer rates compared to the equilibrium boundary-layer theory of Fay 
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Che | — 
tput types of gages, considering the scatter in these types of 10 | : _ . 
at : - 
measurements. * 
e heat- é bs et 
sult lhe theoretical curve on Fig. 12 represents the result 
, | of the equilibrium boundary-layer calculations of Fay |e? | 
eee : 6 
| 
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and Riddell. The equilibrium solution for a Lewis 


































TAVe Is : : le ae a x 
aneous | Number of 1.4 and a Prandtl Number of 0.71 gives a 
¢ OUS | . . ° ° | 
value reasonable correlation to the experimental data within = 
rij the scatter of the data and the uncertainties of the - 7 
) as vite ar . 
lone theory. The uncertainties in the calculated heat => | 
transfer are the uncertainties in the transport proper S| 
é' pee piles Gor a uw | 
+ the ties and the properties in the inviscid flow. The latter Q | 
¢ , . ‘ one | 
“tea depend on the assumptions of thermodynamic equilib- = | 
hy per rium at the stagnation point and a Sutherland de- = DB 
er- | 
) and pendence of viscosity on temperature. 4G 
J, i a us | 
ue rhe consequences of this data can most easily be = 
> from ‘ ; ‘ J 
1 tem. | seen from the nondimensional form of the laminar | 
Leathe | stagnation point heat-transfer equations. The trend | 
/ 
; the heat-transfer parameter. Nu/ Re. wi - L | | 
gravi- of the heat transfer parameter, .\ u V Re, with flight Ol > : : —, 
with velocity, characterized by the density ratio across the eseecee seein, aieitiiiie tiie RM 
: from boundary layer, (pu);/(pu)., Shows a decrease of a factor PH R, HK, 
to be of two in the absolute value of the heat-transfer rate Fic. 13. Effect of the variation of the fluid properties, (pu 
a | from a rate calculated from the low speed value of through the equilibrium boundary layer, on the heat-transfer 
’ " / - 2) ~ a parameter, Nu/~V Re, for Lewis Number of 1.4 and Prandtl 
Nu, V Re from Sibulkin,'® 0.67 at (pu)s/(pH)o = 1.0, Number of 0.71. (Shock tube at initial pressure of 1 cm. of Hg 
as shown in Fig. 13) at the satellite velocity. of air.) 
The Mach Number independence of the stagnation 
point simulation previously described has been con- 
firmed by experiments conducted by the Cornell Aero- 
10% ‘ ‘ yl5 4 ‘ r are co ‘ke > 
nautical Laboratory in a hypersonic shock tunnel. SHOCK MACH NO. -M, 
In this type of facility a gain in flow Mach Number is 2 3 4 5 6 7 8 9 
achieved by expanding air, which has been heated by 19 -—+—F | I ] 
an incident shock wave, in a large area ratio nozzle. |__| it 4 J | | 
= Scat ee ee ae — = as Pc 6 -— DATA REDUCED TO CONDITIONS EQUIVALENT 
By reducing the static tempe rature and speed of sound, — TO STRAIGHT SHOCK TUBE AT P, = Icm Hg | | 
a flow Mach Number of approximately five is achieved. 4 | (APPROXIMATELY 75,000 FEET SIMULATED | ___o$ 
Fig. 14 shows this data, reduced to the same conditions | _ ALTITUDE) AND NOSE RADIUS =.250 INCHES 
and seale as the straight shock tube data, correlated 4 a 
on 2-4. | | Oa 
to the experiments reported here. These measurements, Pp A 
carried out by the Cornell Aeronautical Laboratory 
% - . . . Sc 
under subcontract from this program, were designed 1 lO — : 
to support the present data by establishing the effect ti 08 T 1 
of changing the flow Mach Number from three to five. S 08) 
~ ° » a b + 
Even though this data has been taken at a relatively a 
pel . Fe O4F + 4 
low stagnation temperature only, the measurements “i 
: ° hs ° <a T ; 
are in close agreement with the straight shock tube ry 
data. - 02 Diicenel | | | 
ee ae ee ae ny ee ae Zz THEORETICAL CALCULATIONS 
le range ol altitudes simulated Dy lese experi- S FOR EQUILIBRIUM BOUNDARY 
ments can be extended somewhat if one takes a closer = oj |— _| LAYER L=1.4 
look at the chemical kinetics of the flow over both scale © ost ‘ + | | t ] 
and full size bodies. A simulation based on reproduc Z | 
and full size bodies. A simulation based on reproduc- z 06 | | 
tion of the heat-transfer parameters, Nusselt Number = T 7 t 7 | 
and Reynolds Number, results in smaller actual or | T | 
boundary-layer thicknesses for the small scale model. ] 
The chemical kinetics of the boundary-layer problem —- 02 oa © | | | | | 
a Fe ok de Sh de o PRESENT STRAIGHT SHOCK TUBE DATA | 
1.e., the similarity of the chemical changes in _ the (CALORIMETER GAGE) 
boundary layer—are not simulated completely unless - & CAL. NOZZLE DATA 
the boundary layer is completely in equilibrium  eoeneervesanees 
. FLIGHT VELOCITY x 10° (FEET PER SECOND) 
* More attention has been paid to data from the calorimeter — 14 coo of —" Mach Number on _ stagnation 
= EY ere 1a : i cae einen call ARES point heat transfer on a hemisphere-cylinder body CAL data 
ail instruments because at the high heat tran fer rate f these ex from hypersonic shock tunnel with flow at Mach 4.5 to 5.0 
periments the large temperature rise of the thin gages introduces Present straight shock tube data has a variation of flow Mach 


additional corrections and inaccuracies Number from approximately 2.3 to over 3.0 
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Fic. 15. Schlieren photograph of flow over a flat plate in a 
1'/. in. inside diameter shock tube. The Mach angle should 
be 24.9° if the gas were in equilibrium, 26.2° if the vibrational 
relaxation were very slow, and 31.6° if the gas were completely 
frozen (y+ = 1.4 Measured angles are 25.3 


throughout or completely ‘“‘frozen”’ throughout for both 
scale and full-size bodies. By ‘“‘frozen’’ it is meant 
that the transit time of an atom or molecule through 
the boundary layer to the body surface is very short 
compared to the times required for recombination of 
the atoms. This means that the concentration of the 
chemical species in the inviscid flow at the edge of the 
boundary layer is frozen in the flow and remains con- 
stant throughout the boundary layer. In the inter- 
mediate range, one must consider the simulation of the 
stagnation enthalpy and the recombination rate pa- 
rameter, as used in the solution of the nonequilibrium 
boundary layer with finite recombination rates.’ The 
implication of this simulation is that the mechanism 
of heat transfer by conduction and atom diifusion is 
correctly duplicated. The recombination rate param- 
eter is proportional to the product of the body radius 
and the square of the density and therefore experiments 
performed on small scale models at a given density 
altitude correspond to much higher simulated altitudes 
from the flow chemistry simulation point of view. This 
consideration, when applied to the present experiments, 
0.250-in. nose radius, compared to a body one hundred 
times this size, extends the significance of this data to 
altitudes of almost 200,000 ft. 

Knowledge of the conditions in the shock tube, in 
which the measurements have been made, depends on 
accurate measurement of the velocity of the shock 
wave as it moves down the tube, and on the assump- 
tion that the flow is in equilibrium behind the moving 
The 
velocity measurement brings up the question of atten- 
uation of the strength of the incident shock wave as it 
This question is further com- 


shock wave as well as behind the standing shock. 


moves down the tubes. 
plicated for strong shocks by the properties of equilib- 
rium air, which must be substituted for ideal gas rela- 
Several investigators have devoted much time 
However, 


tions. 
and effort to a solution of this problem.'®!7 
not much work has been extended to the strong shock 
case. 

It is possible to circumvent this question in many 
cases, by making estimates of a tolerable degree of 
attenuation within the desired experimental accuracy. 
Testing time can then by sacrificed by changing the 
shock tube length to obtain the proper compromise. 


AU TICAL 


SCIENCES FEBRUARY, 1958 
Calculations to determine the effect of attenuation 
on the heat transfer measured at a given point in the 
tube have been made, using equilibrium properties for 
the air. This analysis shows that a very small correc- 
tion can be made to the Mach Number, giving an ef.- 
fective Mach Number that would have caused the con- 
the 
general, extrapolating the shock-wave velocity to the 


ditions in which measurements were made. In 
point where the measurements are made is an accurate 
correction. 

Despite these considerations, it has been convenient 
to hold the attenuation for all the data presented in 
this paper to less than 10 per cent of the initial shock 
Mach Number. 
sure where the wall boundary layer is turbulent with 


A few runs made at high initial pres- 


shock attenuations up to 20 per cent have been toler- 
ated. Under this restriction, corrections to the extrap- 
olated shock Mach Number have been calculated to 
be negligible. 

grossly attenuating shock 


Measurements behind 


waves verify this procedure. Both the static pressure 
and heat-transfer rate have been measured to be con- 
stant with time at the value predicted from the extrap- 
olated Mach Number. 

In the shock tube, the question of the thermochemical 
state of the test gas is subdivided into (1) whether this 
gas is in equilibrium behind the incident shock as well 
as behind the standing shock wave at the stagnation 
point of blunt bodies, and (2) whether air starting in 
chemical equilibrium at the edge of the boundary layer 
remains so or is frozen at some composition as it flows 
through the boundary layer. 
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moving shock, made from various assumptions about the thermo- 
chemical state of the air. Note that, in the equilibrium flow 
Mach Number estimate, the vibrational relaxation process across 
the disturbance has been considered 
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STAGNATION 

The assumption that air in a shock tube comes to 
equilibrium very quickly behind the moving shock 
wave is probably very good in the light of available ex- 
perimental data.*” Measurements reported in refer- 
ence 20 of the Mach angles in the flow over a flat plate, 
such as Fig. 15 (a technique similar to the one used 
to obtain pressure distributions), indicate that the flow 
js steady at the Mach Number predicted from equilib- 
Fig. 16 shows the results from 


It can be seen that the experi- 


rium considerations. 
these measurements. 
mental results agree well with the Mach Number cal- 
culated from equilibrium considerations. Any other 
issumption regarding the state of the gas gave dis- 
cernible deviations from the measured angles. 

In addition, a number of other independent experi- 
ments performed at low initial pressure, where the re- 
laxation time should be the longest (including heat- 
transfer measurements), can all be explained only from 
equilibrium assumptions. 

{t can be calculated that, if the flow were not to 
equilibrate across the moving shock wave nor at the 
stagnation point, the effect on the heat-transfer meas- 
urement would be, at the maximum, a decrease of 20 
per cent at the extreme of the range covered by these 
experiments. At less severe conditions of dissociation, 
this correction would be much smaller. 

The question of the chemistry of the air as it flows 
through the boundary layer, starting in equilibrium at 
This point has 
been thoroughly discussed by Fay and Riddel!.* In 
the shock tube, the total heat transfer, whether the 
boundary layer is in equilibrium or is frozen at the edge 
within 


its edge, has received much attention. 


conditions, remains almost unchanged, that is, 
J per cent. For a typical model size, and based on the 
best recombination rate estimates, the large reduction 
of heat transfer that can be produced by a noncatalytic 
wall* in combination with the frozen boundary layer 
can be realized only at initial pressures lower than those 


covered by these experiments. 


CONCLUSIONS 


It has been shown that the environment encountered 
in hypersonic flight through the earth’s atmosphere 


can be created in the laboratory in shock tubes. At 
the stagnation point, these conditions can be used for 
direct measurements of phenomena involving the 


boundary layer, such as aerodynamic heat transfer, 
on scale models simulating larger hypersonic vehicles. 
These conditions have been used to make measure- 
ments of the stagnation point heat-transfer rates on 
hemisphere-cylinder models through a velocity range 
in which the dissociation of the air molecules at the 
Stagnation point goes from zero to over 50 per cent, 
corresponding to satellite velocities. The 
ments are compared to the boundary-layer theory 
developed by Fay and Riddell,® which accounts for the 
large variation of the fluid properties and for heat 


measure- 


* For a discussion of the effects of catalytic and noncatalytic 
‘ay and Riddell. 


body surfaces, see F 


POINT HEAT-TRANSFER 
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transferred by diffusion of atoms in the boundary 
layer. The data and theory show reasonably close 
agreement. The instrumentation and techniques de 


veloped for these measurements are briefly described. 
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Analvsis of Elliptical Rings for Monocoque 
Fuselages 


; STANLEY 


KA U FMA N * 


Aircraft Armaments, Ine. 


SUMMARY 


tae 
prismatic ring subjected independently or collectively to any 


Wise! solved for the external and internal forces of a 


combination of three generalized loading conditions parallel to 


the plane of the ring—a radial force, a tangential force, and an 
applied moment anywhere on the ring 

This paper is a natural extension of bis work as it considers the 
more general shape, the prismatic elliptical ring subjected to the 
following three loading conditions parallel to the plane of the 
first, second, a force 


ring: a force parallel to the major axis; 


parallel to the m.nor axis; and third, an applied moment any- 


where on the ring 


INTRODUCTION 


Fypeseep THERE IS much information available in 
the literature on proposed methods of solving for 
internal forces in a ring of generalized shape*:’ under 
generalized loading, exact solutions have by necessity 
been limited to the simpler geometric shapes such as the 
circle.':'~® For certain specific loading conditions formal 
solutions have also been obtained for the ellipse.’ 

The analysis is subjected to the following assump- 
tions: (1) centroidal axis of the ring and the line of 


shear flow are coincident, (2) strain energy due to 
shear and thrust are negligible compared to that of 
bending, and (3) resisting shear flow is distributed ac- 


cording to the elementary (I’7Q//) beam theory. 
ANALYSIS 
The equation of the ellipse shown in Fig. | is 
y = K'V a2 — x 


Let us now introduce a change in variables by letting 


x=asnr, y = K’acosr 
ds V dx? dy" = av 1 — K? sin* r dr 
where K? + K’”? | 
° ad 
x, =asinz, y, = K’acosz 


The statical moment about the y-axis, Q,, is 
VY, = xds = a[A — A(K,1r)] 


where 


Received April 16, 1957. 
1957. 


* Design Engineer 


Revised and received September 18, 


A(K, r) = | sing V1 — K*sin°¢dd@ 
0 


| _— 1 — K? 
ol — cosr V1 — K®? sin’r) + OK 


(A + 1 
in | = 
(K cosr +4 V1- K 


sin? r) | 
for0 << k <1, A being the complete integral between 
y¥ = Qandr = 7/2. 


The moment of inertia about the y-axis, /,,, 1s 


in = 4 } x*ds 1a°G 
JZ VU 
where> 
G(K, r) = sin? @ V | K* sin? od @ 
JV 
. 
ani 4 K? sin’r sin ¢ cos r + 
(2h? — 1] K’2 
= E(K, 17) +4 ~ F(K, 
3K 3A 


and F(A, r 


the second and first kinds, respectively 


for 0 < k 1, where /:(K, 7) 


integrals of 


are elliptic 














KK, r) = | V1 — Risin? od¢, PUK, r) = 
JV 
[ dd 
woes | K? sin® ¢ 
G being the complete integral between 7 0 and 
w/2. 
The shear stress g can now be calculated as 
g Pizo./1 V1 if (vdx — xdy) | 
J0 
P'4a;+ [A A(K,r)] G + 2coss 7; 
Mt ; is 
—_—_— — 
> M 
N 
¥, 
f x 
a \ 
Ka 
: i ~ ed 
we a — 
Fic. 1. Applied and internal forces, force parallel to major axis 
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Expressions can now be written for bending moment! 
for quadrants I, II, III, and IV in terms of the loading, 
resistant shear flow, and the three redundants .J/o, 7/), 
and V>. One can form, in the usual manner, three linear 
independent equations from the following integrals by 
integrating over the entire ellipse: 

Ol Oo. 0.1/ 


"M ds = 0, "M ds = 0, "M ds = 
J OM, J oll, JS OV, 


After solving for 1/), J), and Vo, equations are ob- 
tained for bending moment (.)/), thrust (7°), and shear 
(V) in terms of the known parameters.’ A_ similar 
procedure can be followed for the other two loading 
conditions and coefficients then calculated for 17, 7 and 
V. The accompanying graphs represent these co- 


efficients. 


CONCLUSIONS 


In the usual airframe construction the centroidal axis 
of the ring does not exactly coincide with the shear flow 
in the skin, attached to the outside flange of the ring. 
It is good practice to choose the greater skin dimensions 
rather than the dimensions of the centroidal axis of the 
ring to be on the conservative side. 

One can resolve any load parallel to the plane of the 
ring into one or more of the three generalized loading 
conditions considered in this paper. Furthermore, one 


. ; ' 
can superimpose the shear flows, bending moments 
shears, and thrusts resulting from any combination oj 
these loading conditions, as the theorem of superpositio: 
for small deflections is valid here. 
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The Effects of Afterbody Length and Mach 
Number on the Normal Force and Center 
of Pressure of Conical and Ogival 
Nose Bodies 


WILLIAM E. BUFORD* 
Ballistic Research Laboratories, Aberdeen Proving Ground 


SUMMARY 


Measurements have been made of the normal force and center 
f pressure of bodies of revolution at supersonic Mach Numbers 
Conical and ogival nose bodies, with cylindrical 


11 calibers, 


from 1.36 to 4.0 
ifterbody lengths varying systematically from 0 to 
The purpose of this study was to describe in a 
of the body aerodynamic 


were studied 


the variation 


systematic manner 
characteristics with afterbody length and Mach Number. These 
jata are in reasonable agreement with the theoretical values 


predicted by Van Dyke’s improved treatment of the first-order 


ixial and cross flows 


SYMBOLS 


dynamic pressure = yP,,.M? 

re = test section pressure 
= ratio of specific heats = 1.405 

VU = Mach Number 
= angle of attack 
Cy = N/qS 
normal force coefficient slope per degree 
normal force 
maximum cross-sectional area of body = xd?/4 
= maximum body diameter 
= body length 
fineness ratio 


C.P = center of pressure 


INTRODUCTION 


yhoo EFFECTS of fineness ratio and Mach Number on 
the normal force and center of pressure of bodies 
of revolution are of considerable interest in missile 
design. Many current missile configurations employ 
long bodies in the design, and in many cases the con- 
tribution of the body to the total normal force is equal 
to or greater than that contributed by the lifting sur- 


The forward center of pressure of the body 


laces. 
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makes it difficult to stabilize such missiles throughout 
the required Mach Number range. 

The art of missile design has been handicapped by a 
lack of systematically acquired information describing 
the aerodynamic variations of bodies with fineness 
ratio, Mach Number, and differing nose configuration. 
The existing experimental data are usually limited to 
describing the effects of one particular variable, or a 
description of a particular missile configuration. Ex- 
tensive theoretical investigations have been made in 
attempts to predict the desired information accurately. 
However, most of the theoretical solutions have suffered 
from inherent defects from the standpoint of accuracy, 
and from the time and effort required to perform the 
calculations. A constructive approach to these limita- 
tions was made by Van Dyke! in a paper in which 
better accuracy was indicated by an improved treat- 
ment of the first-order solution. In the present paper a 
description is presented of the experimental normal 
force and center of pressure variations with several of 
the physical body parameters and with Mach Number. 
In addition to their direct practical value in design, the 
results may be utilized for evaluating the theoretical 


solutions. 


MODELS 


In preparing a systematic program of this type, it is 
necessary to decide which physical parameters are 
most important aerodynamically, and, therefore, which 
revolution the 


ones are to be varied. For bodies of 


parameters we have chosen are (1) nose type (ogive or 
cone) and (2) afterbody length. The parameters were 
chosen and varied in a manner so that the results are 
readily adaptable for use in missile design practice. 
Two basic nose shapes were chosen—a 10° half angle 
cone and a tagent ogive of the same length of the cone. 
The nose was 2.84 calibers long and the ogive had a 
calibers. the models 


and a photograph of the models and 


radius of 8.29 Drawings of 
appear in Fig. 1, 
balance can be seen in Fig. 2. Incremental lengths of 
evlindrical afterbody were added up to a maximum 
fineness ratio of 13.84. 

The region of the nose-body juncture is important 
because of normal force carry-over from the nose to the 


body. Therefore in this region incremental lengths of 
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TABLE 1 
f 8.29R Mach Test Reynolds Number Dynami 
20° + o Number Based on Body Diameter Pressure (psi 

1.36 0.96 & 10° 7.49 
— | 1.72 0.65 X 10° 7.87 

,. OZ 0.25 XK 108, 0.63 10° 2 61,6 
i 2.84 _ 3 55 O.8 X 1080.46 * 10° 1.67,4.17 
Bes: =a 1.0 0.13 & 10%,0.49 & 10! 1.03, 4.2 
- 3.34 - 
— 3.59 - 
- 384 - The angles of attack were c cted for : . 
w 430 “d le angles ol attack were corrected for a strut de 
— 484 - flection due to aerodynamic load and for flow deviatioy 
— 584 - ° i ‘ . = 
Lo 684 ral in the tunnel. In computing center of pressure in th 
' po a small angle range, slopes of the normal force an 
- 11.84 ~ moment curves were used. For the small angles cor 
= 13.84 - 


ALL DIMENSIONS IN CALIBERS 
D=1 CALIBER = 1.490 INCHES 


Fic. 1. Drawing of models. 


0.25 calibers were used. As the body lengthened the 
increments were increased to 2 Bers, These vari- 
ations are presented in Fig. 2. At Mach Number 1.36, 
reflected shock wave interference limited the fineness 


ratio to 7. 
APPARATUS 


The Aberdeen supersonic wind tunnel is a variable 
density, continuously operating type. The maximum 
Mach Number deviation in the test section is +0.015. 
A photograph of the wind-tunnel test section appears in 
Fig. 3. 

The measurements were made with each of two 
strain-gage balances. The balances are shown in Figs. 
4 and 5. The balance shown in Fig. 4 combines a 
normal force beam and moment beam arrangement 
wherein the moment beam is integral with the strut. 
The balance shown in Fig. 5 measures the normal force 
in two parts and may also be used to measure the drag, 


Test CONDITIONS 


Tests were conducted at Mach Numbers 1.36, 1.72, 
3.02, 3.55, and 4. The Reynolds Number was varied 
by changing the tunnel supply pressure. The test 
conditions are summarized in Table 1. The Reynolds 
Number is based on the maximum body diameter (1.49 
in.). In some cases data were obtained up to an angle 
of attack of 16° (see reference 2). However, for the 
purpose of this paper only the data at small angles 
(+35°) are used. Photographs of the flow appear in 


Figs. 6 and 7. 
ACCURACY OF RESULTS 


From a consideration of calibration data, flow uni- 
formity, and data repeatability, it is believed that the 
measurements are accurate to within +0.004 in terms 
of Cy and 0.15 caliber in center of pressure. The nor- 
mal force coefficient slope is believed to be accurate to 
within +2 per cent. A measure of the validity of the 
results is obtained from the close agreement of the cone 
alone measurements with the results of the Stone-Kopal 
exact solution (see Fig. 8). 


sidered, the relationship of Cy versus @ was linear an 
dcy ‘da was obtained by graphical means. 


RESULTS 


(a) Fineness Ratio Effects on Norimal Force 


The variation of the normal force coefficient slope 
with fineness ratio and Mach Number is summarized 
for the ogival and conical nose bodies in Figs. 9 and 1( 
respectively. 
interesting characteristics. 
increase in dcy da as afterbody is added to the nose 
in general, after the afterbody attains a certain length 
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° 
little or no increase results from further extensions. 
For the conical nose, the increase in dcy,/da with addi- 
tion of afterbody length is constant at all Mach Num- 
bers for afterbody lengths up to 1 caliber. At this 
length, deviations occur for different Mach Numbers. 
Also, dcy/da becomes constant at some length de- 
pendent upon the Mach Number. 
configuration exhibits the increase in dcy da 
increased afterbody, and dcy/da@ attains a maximum at 


The ogival nose 
with 


a fineness ratio of about 7 or less for all Mach Numbers. 
However, a Mach immediately 
apparent as the body is lengthened. 


Number effect is 
These data seem 
to positively establish the existence of a potential nor- 
mal force carry-over onto the afterbody; the magnitude 
of carry-over and the length of afterbody affected are a 
function of nose shape and Mach Number. This 
normal force is in contrast with the force developed by 
viscous effects. That is, a normal force is developed 
on the cylindrical afterbody by the potential flow. 


As the afterbody continues to be lengthened, one 
would expect a progressive reduction of the normal 
force carry-over, and this effect seems apparent in most 
of these results. However, there is a possibility of a 
gradual progressive increase in dcy da because of in- 
creased effective base area due to boundary-layer dis- 
placement thickness. This might explain the progres- 
sive small increase in dcy da for both the cone and ogive 
nose bodies at Mach Numbers 3.55 and 4. However, 
this possible effect seems to require more careful study. 
Because of the lower pressure on the upper surface of 
bodies, there may be a tendency for the boundary-layer 
This 


distortion may tend to cancel expected effects of an 


fluid to accumulate on the body's upper surface. 


effective larger base area due to increased boundary- 
layer thickness. In connection with this study, tests 
were made in which, on very long bodies, increased 
displacement thickness was caused by a nose transition 
strip. No measurable differences in dcy da occurred 
with or without the transition strip. 


(b) Mach Number Effects on Normal Force 


In regard to Mach Number effects, these data sub- 
stantiate previous experimental results which have 
shown an increase in dCy/da with increasing Mach 
Number up to a certain limit. This is characteristic of 


pointed nose bodies. 


Any Reynolds Number effect on dcy,da, associated 
with the increasing length, would be included in the 
effect of Mach 
Numbers the Reynolds Number was varied and the 
effects on dcy/da were undetectable. 


fineness ratio. However, for two 


There appears to be a Mach Number where a maxi- 
mum dcy ‘da is reached and further increase tends to 
reduce dcy da. This phenomenon may be distinctly 
observed in the ogival nose data where a maximum is 
reached at 3.02 and the values at 3.55 and 4 are less. 
Also for the cone, dcy/da at Mach Number 4 is smaller 
than it is at Mach Number 3.55. 
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(c) Center of Pressure 


The increase in dcy da, whether caused by increased 
fineness ratio or Mach Number, is associated with a 
rearward movement of the center of pressure. This 


may be observed in Figs. 11 and 12. 


This, of course 
is expected since the increased normal force is con 
tributed by the afterbody. In most cases dcy da and 
center of pressure position become constant at the same 
fineness ratio. One exception is the data for the ogival 
nose body at Mach Number 4; here the center of pres- 
sure location is constant above a length of 8 calibers 
whereas dcy da continues to increase. In actual mis- 
sile practice this rearward movement of the center of 
pressure with increased body length may not be as help- 
ful as it appears since associated with it is the increase 


in the body normal force contribution. 
EVALUATION OF THEORY 


(a) Improved First-Order Theory 


Van Dyke, in a recent paper,' made several important 
contributions, among which are an improved treatment 
of first-order theory, a second-order solution to the 
exact perturbation equations for use in the zero lift 
case, and the suggestion of a hybrid theory combining 
first-order and second-order solutions for the lifting 
case. He notes that in usual first-order theory, several 
approximations are made on the basis that the mathe- 
matical order of the terms neglected is no greater than 
the terms previously dropped in linearizing the equa- 
tions of motion. He demonstrates that the physical 
magnitude of error may be many times larger and that 
while in many cases linearizing the pressure relationship 
agreement is 
obtained with the Dyke's 
hybrid theory has been used to predict the pressure and 


proves disastrous, greatly improved 


exact expression. Van 
normal force on several particular body shapes, and 
quite good agreement has been obtained. It is of 
interest to determine whether the theory has merit 
when applied to more general configurations. 


(b) Solutions Used in Comparison 


The results of the present study provide a criterion 
for an evaluation of Van Dyke’s suggestions. This is so 


because we are concerned both with the flow at the nose 


and also the carry-over onto long afterbodies. Theo- 
retical calculations were performed for the lifting case 
of the variable length bodies with the ogive nose. Two 


differing calculations were used: (1) a first-order 
solution wherein the exact pressure relation and the 
exact tangency condition were applied, and (2) a 
calculation wherein the numerical method of character 
istics was used for the axial flow and a first-order solu- 
tion was used for the cross flow. This latter differs 
from Van Dyke’s hybrid theory to the extent that the 
exact solution rather than the second-order solution 
was used for the axial flow. In all cases, the theo- 
retical values are presented in terms of the normal force 


coefficient, Cy. 
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(c) Evaluation of Improved First-Order Theory 


Evaluation of the improved first-order solution may 
the 


be made from Figs. 13 and 14. In the case of 
normal force coefficient slope, dcy da, the values from 
the improved first-order solution agree within about 
10 per cent of Both the 


observed increase in dcy da@ with increased length and 


the experimental data. 


the maximum value at some particular length are 
predicted. However, at all Mach Numbers, the theory 
predicts a decrease in dcy da when the fineness ratio is 
increased sufficiently. It is conceivable that the theory 
is closely representative of the potential flow, and that 
the constancy of the experimental results above a cer- 
tain length may be attributed to viscous effects—.e., 
boundary-layer displacement thickness. It is also 
noteworthy that, in the ogive nose case, the theory 
predicts the observed decrease in dcy da at Mach 
Numbers greater than 3. 

The predicted centers of pressure are not as good 
quantitatively as are the values of dcy/da; however, 
the qualitative trends of the observed results are 
substantiated. Associated with the predicted decrease 
in dcy da is the forward movement of the computed 


center of pressure above a certain length of afterbody. 


(d) Evaluation of a Hybrid Theory 


The improved first-order and the modified hybrid 
solutions may be compared in Figs. 15 and 16. It is 
apparent that the differences between the first-order 
and the hybrid solutions are of a minor nature. In 
fact, the first-order agreement with the present experi- 
ments is better than that of the hybrid solution, 


although this is probably fortuitous. The predicted 
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much the same, and it that the 


increased effort in using the hybrid solution is not 


trends are appears 


warranted. 


CONCLUSIONS 


In the small-angle, potential-flow range, the normal 
force coefficient slope, dcy/da, increases to some maxi- 
mum essentially constant value as afterbody is added 
to a basic pointed nose. The magnitude of the increase 
and the afterbody length at which no further increase 
occurs are a function of Mach Number and nose shape. 

The normal force coefficient slope increases with 
Mach Number until a maximum is reached. Further 


increases in Mach Number result in a decrease in 
dcy/da. 

The center of pressure moves rearward as afterbody 
is added until a maximum constant value is attained 
and also moves rearward as the Mach Number is in- 
creased. 

The predictions of improved first-order theory are in 
fair quantitative agreement and good qualitative agree- 
ment with the experimental results. 

It is indicated that a hybrid solution 
has little, if any, superiority 


exact axial 
flow, first-order cross flow 


in comparison with improved first-order theory. 
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On the Numerical Calculation of Detached 
Bow Shock Waves in Hypersonic Flow’ 


P. R. GARABEDIAN*® anp H. M. LIEBERSTEIN** 
Stanford Unwersity and the Ramo- Wooldridge Corporation 


SUMMARY 


\ method is described for calculating examples of hypersonic 
flow with a detached bow shock wave past a bluff axially sym 
metric body 
alysis is based on a Cauchy problem for the stream function 
in the subsonic region, where the motion is governed by a pat 
tial differential equation of elliptic type. Through analytic con 
tinuation into the complex domain, the Cauchy problem is re 
formulated in such a manner that it becomes properly set in the 
subsonic region. This leads to a stable scheme for computation 
of the flow by 
stream Mach Number 5.8 are presented in which the flow is de 


finite differences. Numerical examples at free 


termined throughout the subsonic region, and, in particular, the 


detachment distance, the location of the sonic line, and the pres 
sure distribution along the body are calculated 
in excellent agreement with experimental data obtained at the 


California Institute of Technology. 
SYMBOLS 


cylindrical coordinates 


Bt = r,s components of velocity 
y = stream function 
=. pressure 
density 
free-stream Mach Number 
= speed in front of shock 
= speed behind shock 
= ratio of specific heats in front of shock 
= effective ratio of specific heats behind shock 
HH = enthalpy 
S = entropy 
ue = speed of sound 
k = an energy loss across shock 
1 = function of y equal to p/p’ 
1 = temperature 
H = stagnation enthalpy 
= functions of r? describing shock 
= parameters in equation of shock 
c, d = coefficients in partial differential equation 
for y 
P,O = dy/dr, dy /dz 
a, B = characteristic coordinates 
x, 3 = coordinates related to a, 6 
E, 7 = real and imaginary parts of 45 
F = arbitrary function of y 
E,, Eo, E;, Ey = parameters in definition of F 
Ax, An increments in x, 7 
= Avyv/An, ratio of mesh sizes 
o*,a = slopes of the characteristics 
V,Vi, 2, 21,4 = auxiliary functions of ¥ 
B = matrix of coefficients 
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The form of the shock wave is assumed, and the 


These results are 


Ramo-Wooldridge 
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R = vector of unknown functions 

\ = number of mesh points on initial line 

h = distance between body and shock along 
aXIs 

A = curvature of body 

K = curvature of shock 


INTRODUCTION 


I THIS PAPER we analyze the steady supersonic flow 
of an inviscid compressible fluid past a blunt body 
In such a flow, a detached bow shock 
We shall be inter 


ested in studying the properties of the flow in the region 


of revolution. 
wave appears in front of the body. 


between this shock wave and the nose of the body, and, 
in particular, our attention will be focused on the be 
havior of the motion near the stagnation point located 
at the intersection of the nose with the axis of sym 
metry. Indeed, we shall confine our discussion to the 
subsonic region behind the vertex of the detached shock 
wave. Thus it suffices to visualize a model where the 
body extends smoothly to infinity behind the bow wave 
and where no wake or further shocks appear. 

In recent years, a great deal of both theoretical and 
experimental work has been done on this problem (see 
references 1, 3, 5, 8, 9, 10, 11). It has generally 


been assumed that most of the important properties 


and 


of the bow shock wave are of a local character, and the 
above simplified model has proved to be quite successful 
as a point of departure. The most serious difficulties 
encountered in attempting to determine such flows 
theoretically result from the unfortunate fact that for a 
prescribed body the shape and location of the shock 
wave are not known in advance, whence it is actually 
a free boundary problem that must be solved. 

In modern aerodynamics, the Mach Number .\/ ot 
the uniform flow in front of the bow wave must be 
allowed to take on values significantly higher than 


A = &. 


extremely sharp increase in temperature across a strong 


In hypersonic flows of this type, there is an 


shock, and for a real gas it is necessary to take into 
account the One _ possibility 
here is to evaluate the equation of state of the gas by 


effects of dissociation. 
means of tables of the pressure p, enthalpy //, and en 
tropy S as functions of the temperature 7’ and density 
p which have now become available*® for a wide range 
7 and p. However, in our analysis the 
situation simplifies considerably because we restrict 


of values of 


our attention to the subsonic region near the stagnation 
point. Crude arguments show that in this region / 


and p vary by as little as a sixth of their total values, 
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JOUR 
so that only a very small section of the tables is actually 
relevant for the subsonic portion of a specific flow. 
Examination of such a section of the tables establishes 
that, over the small ranges of values of 7 and p that 
we have described, the quantity //p p deviates from an 
appropriate fixed value by less than | per cent. Thus 
we are justijed in using in the subsonic region the em 
pirical equation of state 


IT ly (4 1)] (Pp p) (1) 


where y is a constant to be found from the tables and 
can be interpreted as an edective ratio of specific heats. 

In practice, y always exceeds |, but must be chosen 
quite close to | when the temperature in the subsonic 
region is high. In front of the bow wave, of course, an 
actual ratio of specific heats y; replaces the empirical 
value y in the equation of state (1). It is also desirable 
to provide for the addition of a constant on the right 
in Eq. (1) in case this should result in a better repre 
sentation for the enthalpy behind the shock wave. It 
can easily be shown from the principles of thermody 
namics that this slightly more general form of the equa 
tion of state is equivalent to the assumption that the 
entropy S is a function of p p’ alone. 

Determination of the shape of the bow wave in front 
of a prescribed body involves, as we have observed, the 
solution of a difficult free boundary problem. On the 
other hand, it is possible to find the body generating 
a given shock wave by solving an initial value problem 
for the equations of motion. Such an initial value prob 
lem is not properly set in the usual sense when we assign 
Cauchy data along the are of the shock wave bordering 
the subsonic region in a meridian plane, since in this 
region the flow is described by a stream function which 
satisfies a partial differential equation of elliptic type. 
However, we shall show how this initial value problem 
can be reformulated in the complex domain in such a 
manner that it becomes properly set, and our reformu- 
lation will lead to a stable numerical scheme for calcu- 
lation of the stream function in the subsonic region. 
This procedure enables one to obtain numerical ex 
amples of the flow with a detached shock past a blunt 
body of revolution. 

Such examples yield much significant information 
about the physical properties of the flow. It is pos 
sible to deduce from them a theoretical graph of the 
product xh of the curvature « of the body times the de- 
tachment distance / as a function of the free-stream 
Mach Number J/. A description of the distributions 
of pressure and density along the body can also be 
found, and the position of the sonic line can be deter- 
It is reasonable to expect that all this data is 
to a certain extent insensitive to alterations in the 
shape of the body, and thus the reliable examples 
which our method gives will be of considerable practical 


mined. 


value. 


CAUCHY PROBLEM FOR THE STREAM FUNCTION 


We shall work in a meridian plane with the usual 
coordinates r and z, where 7 is the radius and ¢ the alti- 


TIC 
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tude of a cylindrical coordinate system in space. Ij 
we denote by uw and v the components of velocity in th 
directions of increasing y and 3s, respectively, then the 
conservation of mass implies the existence of a stream 
function WY such that 


y, = rpv, wy 


Vr pu ? 
The conservation of energy can be expressed by means 


of either of the two relations 
HT(p, p) + (97/2) I] 3 
p A(wW)p t 


where // is the function of p and p defined by Eq. (1 
77 uw” + v* is the square of the speed, /7p is a constant 
called the stagnation enthalpy, and A is an arbitrary 
function of Y alone. The equations of conservation of 


momentum can be combined with Eq. (5) to give Eq 


(4), and conversely, they can be combined with Eq. (4 
and (+) to elimi 


to give Eq. (3). Using Eqs. (2), (3), 


nate derivatives of u, v, p and p, we can derive from 


the momentum equations the second-order partial 
differential equation 
[1 — (u?/u?)]v,, — (2uv/u*) p,- + [1 (v? uw”) ly 

pu — r°o*[(FT + 9?) y] [A’(W)/A(y)] (5) 
for the stream function y, where yu” yAp is the 
square of the speed of sound. The conditions (2), (3), 


(4), and (5) comprise our formulation of the equations 
governing the steady axially symmetric motion of a 
and (+4) can 


compressible fluid. Since Eqs. (2), (3), 


in principle be solved for u, v, p, and pin terms of 7, y, 
y,, and y., this system reduces, in effect, to the single 
second-order quasi-linear equation (5) for y, which is 
of course, of the elliptic type in a region of subsonic 
flow. 

Without loss of generality, we can choose a scale 


such that p = p 1 in front of the bow shock wave, 


and we denote by q \Jy, * the speed there. We 
write the equation of the bow wave in the form 
(2 + 20)? = f(r*) = 207 + g(r?) (0 


in the meridian plane, with g(0) = 0. The shock con- 
ditions across this curve which express the laws of con- 
servation of mass and momentum are the usual ones of 
classical fluid mechanics, but, because we want to take 
into account the effects of dissociation in our analysis 
we must write the shock condition expressing conser- 
vation of energy in the generalized form 


[y/(y — 1)] (b/p) + (q?/2) 
ly/(y — 1)] + (qi2/2) + Rk (7) 


where /, p, and g stand for values of the pressure, density 
and speed just behind the shock and where & is a con- 
stant which allows for the change in the formula for 
the equation of the state as we cross the shock. Thus 
k depends on y: and on the additive constant which 
might be included in Eq. (1). 
where no additive constant is inserted in the empirical 
equation of state (1), we find for k the value 


In the simplest case 
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NUMERICAL CALCULATION 


(y — 1)] (S) 
[he constant on the right in Eq. (7) is, of course, merely 
another expression for the stagnation enthalpy JJ». 
Since the flow in front of the bow wave is uniform, 
ind since the stream function y is continuous across the 
shock due to the conservation of mass, we find that 


y= alr? 


2) (9) 


on the curve (6). Straightforward calculations with 


the remaining shock conditions show that along Eq. (6) 


Y: (gir2f'(r?2) (2 + o)LJ [1 -— A V)] (0) 
ind that 
1(y) ty = 1)/( 1)] X 
(Li — gx? VEL — VL?) L?] V (11) 
where 
L= 14+ v2 m)? af a) (12) 


Ly = G2k + [y/(y — DIL + aL 13) 
Aji — + CL qi) 14) 
. (y — 1) (qi? + 2kL) + 2yL 
| i172 (15) 
ya + 7yL + gilgi7A* — Ay = 1)kL | ~ 
We note also that 
Ay) L(y) 
—" 
A(y) ] 
7 = 1L)LAV; - y¥ViVL (91? T L) _ qa? V1 - ; 
= (16) 
Ly i VL(g1" L) 
where 
1" YY 
oe: l)q 
4 HAF CY? = Dae 2 
z ) (la) 
| [gi2A2 — 2(-y2— 1)RL]'/79 


We can now rewrite the second-order partial differ- 


ential equation (5) for y in the form 
+ ¢W, 


and d are known functions 


aw + 2b y +d=0 (1S) 
where the coefficients a, 2, c, 
ofr, p, v,, 
of the relations (2), (3), (4), (11), and (16). 


and W. which are defined explicitly by means 
When the 
curve (6) is given analytically, Eqs. (9), (10), and (18) 
constitute a Cauchy problem for the determination of 
the stream function y. With our choice of the bow 
wave in an appropriate class of analytic curves (6), we 
shall investigate this Cauchy problem in the subsonic 
region, where Eq. (18) is of elliptic type. In particu 
lar, we shall work through in detail the case 


(\/? — 


1) (22 + 


a? = [?(M? — 


1) (19) 


in which the shock is a hyperbola whose asymptotes are 
Mach lines of the approach flow. 
CHARACTERISTIC COORDINATES 


The two families of characteristics of Eq. (18) are 
defined by the quadratic ordinary differential equation 


OF 
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az? — 2b2% + cf 0 20) 

This equation can be factored to yield 
(3 — ott) (8 — oF 0 21 


and for the moment we shall not let it disturb us that 


the roots ¢* and o~ are complex numbers when / 
ac <0. Each factor in Eq. (21) defines a family of 
characteristics of Eq. (1S). We set 

P=y, Q v 22 


and we introduce new independent variables @ and 3, 
called characteristic coordinates, which are defined by 
the property that one of them is constant along each 
curve of one family of characteristics, while the other is 
constant along each characteristic of the remaining 
family. 

In terms of the characteristic coordinates a and 8 
(18) 4 


we can reduce the second-order equation in a 


standard way” * to the system of six first-order partial 


differentia] equations 


_& oO r (), c } 0) 23 
datr,taa?* P, ca, 0 24) 
do-tataca Ps, C Os 0 25 

VY. — Pr. Oz. 0 26) 
W, — P? Oz 0 97) 


for the determination of the five unknown functions 
r,s, P,Q, and y. We do not need to use both of Eqs. 
(26) and (27), and, indeed, it suffices to take any fixed 
linear combination of them, although such a combina 
tion introduces in general a new artificial family of 
characteristics in the (a, 8) plane. Together with the 


initial conditions 


zg = g(r) [f(r?)!/2 + ay), y qr*/2 IS) 

QO [qir2f’(r?) (2 + 2o)L] [1 ifs 29) 
P=q rf'(r?)O (s + 2 30) 

the first-order system (25), (24), (25), (26) is equivalent 
to the Cauchy problem defined by Eqs. (9), (10), and 


(18). Since further transformation of the characteristic 
coordinates is still permissible, we can impose the initial 
conditions along an arbitrary monotonic curve in the 
(a, 8) plane and prescribe r arbitrarily along this curve 

It turns out to be convenient to make the additional 


substitution of variables 


t= (e+ 6)/2, 3 (a 3) 27 1) 
Eqs. (23), (24), and (25), together with the sum of 
Eqs. (26) and (27), transform into the system 
, 7 l ‘ 
Yr, = |07, az ac = oF)" 32) 
. . ; 19) 1/2 
Ss 1c ry Zi (ae = 0*) > 
r. = (d Sy T b P, 7% QO (a¢ bh? i 54) 
i, = [-—d fo — @ F. —h QO, ac — b*)'’* 30) 
vy, = [(Pb + Qc)r,— (Pa + Qb)z,] (ac - b2)!* 36 ) 
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CHARACTERISTICS 
4 
Fic. 1 Geometry of characteristic coordinates in the complex 
domain 
in the (x, y) plane. This canonical system, together 


with the initial conditions (28), (29), and (30), is also 
equivalent to the Cauchy problem (9), (10), (18), and 
it has the advantage that its real solutions generate 
real flows in the subsonic region, where 6” — ac < 0, 
since no complex quantities appear in that case. 
the sake of simplicity, we rewrite Eqs. (32), (33), (34), 
(35), and (36) in the matrix notation 


R, = BR, 


(37) 


where RX stands for the column vector of unknown func- 
tions r, s, P, Q, wy and B is the appropriate five-by-five 
matrix of coefficients. 

We are free to choose the y axis as our initial curve 
in the (x, vy) plane, and we impose there the conditions 
(28), (29), and (30). To these we can add a relation 7 = 
F(y), where F is an arbitrary real analytic function of 
y. This freedom results from the easily verified fact 
that our formulation of the Cauchy problem is invariant 
under any conformal transformation in the (x, y) plane. 
For the applications, it is convenient to use an explicit 
initial relation between ry and y of the form 


r= F(y) = Evy + Eny® + [Eny/(y? + Ey)] (38) 


where /4), Eo, Es, 


and /; are arbitrary real parameters. 


UsE OF COMPLEX VARIABLES 


In our reduction of Eq. (1S) to the canonical form 


yQ> 


(37) we have not allowed the appearance of complex 


quantities to disturb us. This is justified when we 
take into account the analyticity of the functions which 
occur and when we think of the intermediate equations 
in our derivation as identities which must hold in the 
complex domain, according to the basic principles of 
analytic continuation. Thus we should assume from 


the beginning that we have continued all functions 
analytically into the complex domain so that they have 
It is 


then irrelevant whether or not complex numbers are 


become complex functions of complex variables. 


used in our calculations. It is convenient to observe at 
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the conclusion of the argument, however, that, whe; 


b? — ac < OQ, real solutions of the canonical system (37 


in the real (x, y) plane correspond to real subsonic flows 


satisfying Eq. (18) in the physical plane. 

We shall now make a more decisive use of comple, 
variables by showing how the Cauchy problem for Eq 
(18) can be formulated in the complex domain as a prop 
erly set initial value problem for the canonical systen 
(37). For this purpose we keep x real, but we let y = 
£ + i be complex. 
of the independent variables in a space of three dimen 
sions with the rectangular coordinates x, &, and yn. The 
initial conditions (28), (29), (30), and (38) are imposed 
for x = 0, and they become identities in the complex 
variable y over the (&, 7) plane. In particular, they ar 
valid for variable 7 with x = O and & fixed. With this 
in mind, we shall restrict our attention to a plane & = 
const., in which x and 7 are the rectangular coordinates 

Because of the Cauchy-Riemann equations, we car 


write Eq. (37) in the form 


RK, = (6/12)K. 29 


in such a plane. The matrix B 7 is easily seen to have 
the characteristic roots 1, —1 and 0, and the system 
(39) is of hyperbolic type in the (x, 7) plane and has 
the three families of characteristic curves x + 7 = 
const.,.. — 7 = const., and 7 = const. The last family 
appears because we used the sum of Eqs. (26) and (27 
in our derivation of Eq. (37), whereas the earlier tw 
families correspond to the original characteristics of the 
In this 


connection, it is interesting to note that Eq. (51) re- 


second-order partial differential equation (18). 


duces, at least in the special case & = 0, to the rotation 
and magnification 


») _ 


(a+ 8)/2, n= (8 — a) 2 (40 


x = 


and is therefore a real substitution in the (x, 7) plane 


Thus our introduction of 7, rather than y, as the basic 
variable is equivalent to using real characteristic co 
ordinates. (See Fig. 1.) 

What we have established at this stage is that, by 
analytic continuation into the complex domain and 
appropriate choice of the independent variables, we 
can transform the real elliptic equation (18) into the 
hyperbolic system (39) for the determination of five 
complex-valued unknowns. Furthermore, we have 
shown that when 7 = 0, that is, in the (x, &) plane, the 
solution of Eq. (39) will be real. 
the unknowns ¢, z, P, Q, and y will in general turn out 
to be complex because of the appearance of 7 in Eq 
(39) and because of the complex values of y = & + 
which we substitute initially into Eq. (38). 

In order to find an actual subsonic flow in the physical 
plane, we must solve the initial value problem for Eq 
(39) for choices of the parameter £ covering a complete 
interval. Thus with our method we must determine 
the solution in a three-dimensional region, although 
only the portion of this region which intersects the 
plane 7 = How- 
ever, at each specific stage of the procedure, we have 


0 yields physically significant data. 


Thus we can visualize the domai 
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NUMERICAL CALCULATION 
merely to deal with an initial value problem in two 
independent variables. The solution of each individual 
problem of this type yields values of y, P, and Q along a 
curve in the (7, z) plane. This curve assumes a shape 
which depends on our choice of the arbitrary function 
Fin Eq. (38) and which can be prescribed quite gener- 
ally 
The Cauchy problem for Eq. (39) is properly set in 
the sense that the solution depends continuously on the 
initial data, but the same is not true for Eq. (18) in the 
elliptic case b> — ac <0. The explanation for this phe- 
nomenon is that a small change in the initial data for 
Eg. (1S) in the real domain can result in a large change 
in the initial data for Eq. (39), since we obtain the latter 
by analytic continuation of the former into the complex 
Thus we should treat the Cauchy problem 
(1S) only when the data can be continued 


domain. 
for Eq. 
into a sufficiently large section of the complex domain. 
The problem becomes properly set if we adopt a new 
criterion for continuity in the variation of the data 
which states that two sets of data are close only when 
the analytic continuation of their difference is small 
in the complex domain. Hadamard’s classic counter- 
example does not satisfy this criterion. 

An even more precise analysis of the situation can be 
made if we are content with a discussion of what hap- 
pens in the space of the rotated characteristic coordi- 


nates x, &, and 7. Here singularities are observed to 
travel along the characteristics x + » = const. and 
v — n = const. Hence a singularity of the initial data 


at the point with coordinates £ = £, 7 = mo in the plane 
x = O will appear in the real (x, &) plane as a singularity 
of the solution at the two points whose coordinates are 
Thus singularities 


v = no, & = &, andx = —7n, — = &o. 
of the initial data which occur for relatively small values 
of n should be avoided. Also, if we wish to alter slightly 
Eq. (6) of the given bow shock wave so as to approxi- 
mate more accurately the flow past a prescribed ob- 
stacle, we should only make changes which remain small 
after analytic continuation into the portion of the (&, 7) 
plane that projects along the characteristics onto the 
region of the (x, €) plane corresponding to actual flow. 

We can now justify more rigorously our choice of the 
hyperbola (19) as an appropriate shape to take for the 
bow wave. In the case of the hyperbola, indeed, the 
only singularity whose existence is predicted by our 
method would occur at the focus if Eq. (18) were to 
have the same fixed complex characteristics as Laplace's 
equation, since the characteristics through the focus 
would then be tangent to the initial curve in the com- 
plex domain. As matters stand, the singularity should 
occur not too far from the focus, whereas we can expect 
to obtain a flow past a smoothly rounded obstacle which 
lies much closer than this to the hyperbola. 

Our analysis indicates that quantities of physical 
significance such as the detachment distance /; between 
the body and the shock wave depend on properties of 
the flow over the whole subsonic region. Thus we 
could not hope to obtain a reliable description of, for 


example, the dependence of / on the free-stream Mach 


OF 


DETACHED BOW SHOCK WAVES 
Number .\/J by a method which, unlike ours, did not use 
Cauchy data along the complete are of the bow wave 
bordering the subsonic portion of the flow. Our tech 
nique based on the canonical system (37) determines 
the flow throughout the subsonic region, but not be 


” 


4 3), (34), 


yond, since the denominator in Eqs. (2), 
(35), and (36) vanishes on the sonic line. The location 
of the sonic line itself, however, can be found by using 
Eq. (37). In the supersonic region beyond the sonic 
line, the classical procedures for treatment of the initial 
value problem for a hyperbolic equation apply, and 
therefore we shall not concern ourselves in this paper 
with the apparently more routine difliculties which 


present themselves there. 


FINITE DIFFERENCE SCHEME 


We are interested in giving a numerical solution of the 
(30), 
dit 


initial value problem defined by Eqs. (28), (29), 


(38), and (39). We shall outline briefly a finite 
ference procedure which is well suited for this purpose. 
We denote by Av and Ay positive increments in the 


variables x and 7, and for a given fixed value of € we set 


Rn, n = 
Bus a = 


R{(mAx 2), (nAn 2)]\ 


11 
B[(mAx 2), (nAn 2)]9 ves 


where R(x, 7) and B(x, n) stand for the vector and ma 
trix appearing in Eq. (39), thought of now as functions 
of x and y. For our numerical analysis, we approxi 
mate the partial differential equation (39) by the differ 


ence equation 


} ae = R,, La TT @ x 
[te im ae 1) 211 (Res, wit _ a 1) 


})) 


An. Since the characteristic roots of the 
matrix B i are 1, —1, and O, the solution of the initial 
value problem for Eq. (42) is stable if and only if s < 1. 
the 


where s = Ax 


Because central differences are used, truncation 
error is of the order of magnitude of the square of the 
mesh size An. 

In practice, it is only necessary to use values of m and 
n in Eq. (42) whose sum is even. Furthermore, since 
R(x, 0) and B(x, 0) have real elements, the Schwarz 
principle of reflection shows that 


x... , = z.. ls B,,. 1= B, “ (43) 


where R&,,,, and B,,, ; denote the matrices whose ele- 
ments are the complex conjugates of the elements of 
R,,,; and B,,,1. This symmetry property permits us 
to solve the initial value problem for Eq. (42) using 
only mesh points with » > 0, since when n > | we do 
not need any mesh points in the lower half-plane in 
order to evaluate Eq. (42), while for » = O the mesh 
points in the lower half-plane can be eliminated from 


Eq. (42) by means of Eq. (43). Thus our numerical 


scheme will consist first in calculating Ro, , for n = 0, 
9 4,...,2N — 2 from the initial conditions (28), (29), 
(30), and (38), then in calculating Rj, ,, for m = 1, 3, 5, 


...,2N — 3 by means of Eq. (42) with R_4, , replaced 


by (Ro, nut + Ro, n-1)/2 


and s replaced by s 2, and 
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finally in calculating through repeated application of 


Eq. (42) the remaining values of R,,, , such that m + n 
is even and does not exceed 2V — 2. (See Fig. 2.) 
The method which we have now presented for the 
determination of hypersonic flows with a detached shock 
past bluff bodies of revolution is based entirely on the 
sound principles of exact theory. It can hardly be 
considered to be elementary, however, since we have 
been obliged to introduce complex quantities in a 
physically irrelevant space of three dimensions. We pro- 
ceed to demonstrate that the approach is nevertheless 
quite practical by giving explicit numerical examples 
computed in all detail on a high-speed electronic calcu- 


lator. 


DISCUSSION OF NUMERICAL COMPUTATIONS* 


The procedure we have outlined in the previous sec- 
tions has been programed for the ERA 1103 electronic 
calculator at the Ramo-Wooldridge Corporation. The 
program is designed so that the physically significant 
parameters 7, y, and y; can be varied freely as well 
as the geometrical parameters /, &, fy, Fo, A:, 4 and 
the mesh size An. For a specific choice of these param- 
eters, a typical run with the number JN of grid points 
on the initial line equal to 36 takes less than one hour 
of machine time. 

Alterations in / merely produce a change of scale in 
the physical plane, while through variation of & we can 
sweep out curves of data there which cover a complete 
region. The influence of the parameters /,, Hs, £3, and 
kk; on the calculations is more subtle, and we shall 
subject it to a detailed analysis. 


Let us suppose that we have carried out the compu- 


tations for a series of values of ~, with EF, = bk; = 0, Fy 
= Kk, = 1, and with the remaining parameters held 


fixed. We are interested in estimating the changes 
which can be effected by introducing new values of Fi, 
ky, E3, and Fy in Eq. (38). The original computations 
establish a relationship between points in the real 
x, ¥) plane and points in the physical plane. Now an 
alteration in the form of the function F inserted in the 
initial condition (38) can be interpreted as the result 
of performing the conformal mapping Y — iY = 
F(y — 1x) in the (x, y) plane, since our derivation of the 
canonical system (32), (33), (34), (35), (36) is invariant 
under such conformal transformations. It follows 
that we should be able to determine the correspondence 
between points in the transformed (x, y) plane and 
points in the physical plane directly, without additional 
use of the system of partial differential equations. 
Indeed, if the point (X, Y) is the image of the point 
(x, y) under the conformal mapping Y — 1X = F(y — 
ix), and if the original solution of Eq. (39) with £,; = 
FE, = 1, Ey = BE; = 0 gave at (X, Y) the pair of coordi- 
nates (7, z) in the physical plane, then the solution with 


* The authors wish to acknowledge their gratitude to Robert 
Beach, Kenneth Curewitz, Lynn Stoller, and Clyde Coller for 


their assistance in the preparation of this program. 
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the new initial data r = F(y) should yield these same 
coordinates (7, 3) precisely at (x, y). 

Since different choices of F correspond to different 
choices of a conformal transformation in the plane of 
the canonical coordinates (x, y), it is evident that for a 
given value of £ we should be able to obtain a variety of 
curves of data in the physical plane by appropriate 
selection of the parameters /;, 2, -;, and Ky. The 
solutions of Eq. (39) which we arrive at in this manner 
generate widely differing two-dimensional surfaces in 
the four-dimensional complex extension of the (7, z 
plane. In the numerical applications of our method, 
we require an F such that the curves of data which 
sweep out the physical plane as & varies will not only 
cover the body behind our prescribed shock wave, but 
will also reach the sonic line. To achieve this objective, 
we have taken 


F(y) = (Ey + y3/9)/(E + y?) (44) 


corresponding to Eq. (38) with £; = 1/9, A, = 0, E 
SE 9, and Fy = Ek. When F is given by Eq. (44), the 
images of the lines y = const. under the conformal 
transformation Y — iX = F(y — ix) are curves wind- 
ing around a double critical point in the (XY, Y) plane 
aA = Y = (2/3) We choose F so that this 
critical point maps into the sonic point on the shock 
wave in the physical plane, whence our curves of data 
bend around to meet the sonic line. 


The expression 


F(y) = Ey/(E + y?) (45) 


yields another useful choice for / which results from 
Eq. (38) when we set £; = FE. = 0, FE; = Ey = E. Since 
in this case the critical point in the (Y, Y) plane, which 
occursat X = 0,Y=E 
obtained which bend less sharply than those generated 
by Eq. (44). 

So far our discussion of selecting parameters has been 


2, is simple, data curves are 


confined to the problem of obtaining adequate data be 
hind a fixed shock wave. We shall consider next a more 


general shock curve than Eq. (19), which is defined by 
f(r?) = (1, ak lar”), Zo" Lil (46) 


This equation describes a conic section which can be 
either an ellipse or a hyperbola. We note that in the 
case of the ellipse, however, the formulas in this paper 
yield velocities which must be reversed in order to gen- 
erate a properly oriented flow. We wish to establish 
here a qualitative analysis of the dependence of the 
shape of the body behind the shock (46) on the param 
eters /;, /s, and /3. 

Our conclusions in this direction are based on the re 
mark that singularities of the flow can be expected to 
occur near the foci of the shock curve. Indeed, we have 
seen that if the stream function were to satisfy La 
place’s equation instead of Eq. (18), a characteristic 
passing through a focus of Eq. (46) would be tangent to 
this initial curve in the complex domain. Thus, even 
in the case of Eq. (18), we anticipate that the stream 
function will be multiple-valued in the neighborhood of 
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the focus and will have a branch point near by. A 
branch point is exactly the type of singularity exhibited 
by flow atacorner. Hence we should be able to achieve 
corners on the body behind the shock wave (46) by 
placing the foci of this conic section in the positions 
where we desire that the corners should appear. Thus 
to obtain flow past a rounded cone we should choose 
);, /2, and /; so that Eq. (46) is a hyperbola with its focus 
just behind the stagnation point, whereas to obtain 
flow past a flat-nosed body of revolution we should 
take Eq. (46) to be an ellipse whose foci are located off 
the axis of symmetry near the points where the shoulder 
intersects the meridian plane. It is hoped that these 
comments might be helpful to investigators who wish 
to select curved shocks designed to generate a pre- 
scribed family of body shapes. 

In this paper we present numerical examples which 
serve to establish the success of our method in practice. 
We have made our sample runs at free-stream Mach 
Number 5.8 with a ratio of specific heats y = y. = 14 
in order to facilitate comparison with experimental data 
obtained in the hypersonic wind tunnel at the California 
Institute of Technology. The correlation between the 
theoretical and the experimental determinations of de- 
tachment distance and pressure distribution along the 
body is excellent, as is evident from a check of our 
tables and graphs against the results given by Oliver.’: '° 

Our principal calculation has been made for the 
hyperbolic shock (19) with / = M = 5.8, ky = 1/9, 
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Ek. = 0, BE; = 0.5163, EA, = 0.5808, and we have at- 
tempted to tabulate here enough of our data for this 
case to give the reader an accurate description of the 
information obtained. In the actual computations 
complex floating-decimal arithmetic retaining eight 
digits was used. However, in our Tables we have 
rounded off the results so that only as many figures are 
listed as have been shown to be significant by trunca- 
tion error tests based on variation of the mesh size, 
Roughly speaking, we had to take V = 532 grid points 
on the initial line in order to ensure getting four sig 
nificant figures in the data along the body. This con 
vergence deteriorates as we approach the sonic line be- 
cause of the smallness there of the denominator (ac 
b*)’ in Eqs. (32), (35), (34), (35), and (36). Thus our 
results along the sonic line have an accuracy of only 
two or three significant digits. 

In Fig. 3 we present a detailed diagram of the flow just 
referred to. In this Figure the positions of the body 
and of the sonic line are shown, and velocity directions 
are indicated at sequences of points where runs of our 
program have punched out the flow data. The loci 
described by these points have nothing, of course, to do 
with physical properties of the flow, but depend rather 
upon the geometry associated with our method. We 
also include a graph of the pressure distribution along 
the body (see Fig. +) and insert for purposes of com 
parison the corresponding experimental values due to 
Oliver® '° for pressures on a sphere. The agreement is 
good when we take into account the differences in body 
shape between his work and ours. 

Figs. 5 and 6 show two further diagrams of flows 


computed using a more general shock curve of the type 


(46). Fig. 5 corresponds to the values /, .-== 6. 
I; = 1, M = 5.8, y = ym = 1.4 of our principal param 
eters and depicts flow in front of a rounded cone. Here 
all runs were made with £, = E. = 0, FE; = Ey, = 0.9746, 
except for one with £; = Fy, = 0, Fy; = Ey = 1.456. 
Fig. 6 was obtained with /; = 1.095, /2 = 0.0865, /; 

—1, W = 5.8, y = yn = 1.4, and, for the most part, 


with E; = ky= 1, k,2 = E;=0. It represents a flow past 
a body with a quite flat nose. Our purpose is to dem 
onstrate with these Figures how readily our analytical 
technique permits us to calculate hypersonic flows with 
a detached shock past a variety of bodies of the kind 
in common use. Values for the pressure distributions 
along these bodies have been presented directly on the 
diagrams in order to make possible a comparison with 
experimental data. Again, our results agree quite 
well with those compiled at the California Institute of 


Technology hypersonic wind tunnel. 


EXPLANATION OF THE TABLES 


In Tables 1, 2, and 3 we list data for the flow we have 
calculated with 1 = 5.8, y = 7, = 1.4,8 = 0, s = 1, 
1 = 538, E, = 1/9, &, = 0, E, = 05163, EF, = 0.5808. 
Since we have normalized so that the pressure /; and 
density p; in front of the shock are both unity, the 
values of p and p that we present are actually to be 
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) 


PABLE 3 


Data on Sonic Line for Fig. 3 


interpreted as values of pp; and p p;. Similarly, the 
velocity components uw and v as we list them are really 
values of the dimensionless ratios (p:/fi)‘ u and 
(p: pi)’ v. Table 1 represents flow quantities on the 
shock, in the interior of the shock layer, and, when y 
is negative, inside the body. Excerpts from each run 
of the program corresponding to a specific choice of é 
and An are listed under a single heading, and the runs 
are numbered consecutively in the order of increasing 
— With each entry is included the value of the index 
m in our finite difference scheme for which it was com- 
puted. The data in Table 2 is evaluated for points on 
the body by interpolation, and the run from which each 
point was obtained is identified. Table 3 was derived 
by extrapolation in a similar manner and describes the 
flow at points along the sonic line, including one point 
inside the body. 

Fig. 3 is a diagram of the flow described in more de- 
tail in the tables. Velocity directions at a selection of 
points in the interior of the flow are marked, and the 
product of the detachment distance / times the curva- 
ture « at the vertex of the body is listed. In Fig. 4 we 
compare graphically plots of the pressure distribution 
as a function of the angle of inclination of the body for 
this flow and for Oliver’s experimental work with a 
sphere. 

Fig. 5 shows the flow in front of a rounded cone ob- 
tained by running our program with 1/7 = 5.8, y = y1 
l4,4k=O,s=1h=h =08,)5 = 1, £, = E, = 0, 
and with #; = £y = 0.9746 or E; = Ey = 1.456. The 
positions of the body and of the sonic line are shown, 
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together with flow directions at interior points, and we 
have included in the diagram values for the pressur 
along the body. Fig. 6 is a similar diagram for flow 


past a flat-nosed body of revolution and was derived 


by setting 1/ 5.8, ¥ v1 1.4, k u, 5 
i 1.095, Io 0.0865. / 1, and, primarily, / 
Es l.£:= 2 0 in our basic program 
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Effect of Quasi-Steady Air Forces on 


Incompressible Bending-Torsion Flutter’ 


J. DUGUNDJI* 
Massachusetts Institute of Technology 


SUMMARY 


Explicit solutions are obtained for the bending-torsion flutter 
if two-dimensional airfoil in incompressible flow under the 
issummptions that the Theodorsen function, C(k) is set equal to a 
real constant, and the diagonal virtual mass terms are negligible. 
For the case of small bending to torsion frequency ratio, w;/wg, a 
comparison is made of these quasi-steady solutions with an 
earlier empirical expression suggested by Theodorsen and Garrick 
for the nonsteady case, and the effect of the C(k) function is 
indicated. The importance of the c.g. location for these small 
cases is re-emphasized, and the possibility of flutter at zero 


iir speed is indicated. 


DISCUSSION 


ge EFFECT of aerodynamic lags as expressed by the 
Theodorsen function C(k) plays an important role in 
flutter of airfoils in incompressible flow.':? It might be 
of general interest to examine this effect and indicate 
its magnitude for certain specific cases. 

First, the flutter of an airfoil using quasi-steady air 
forces is to be determined (see, for example, Chap. 6 of 
reference 5 It has already been shown that in the 
high supersonic case where a quasi-steady theory (piston 
theory) is valid, rather simple closed form solutions 
result for the two-dimensional bending-torsion case.‘ 
Following a suggestion of Prof. H. Ashley, the incom 
pressible quasi-steady theories are examined in a similar 
fashion to observe their forms and trends. 

The equations of motion for the bending-torsion 
flutter of a two-dimensional airfoil in incompressible 
flow can be written in the Laplace transform, non 


dimensional form, 


(7 ] of DL, V7 
{iu + 1]p? + 2C(p)p + wlan we) kath 4 


(utp? + {1+ 2[(1/2) —a]C(p) ip + 2C(p))ae=0 (1 


[ME gp 2[(1 2) + a]C(p)pyh 4 
([ure a (1/8)]p? + , (1/2 a 
oT Va ) 9 
2{(1/4 — a’ |( (P)yPpr Urea 
2fa + (1 2)|C(p) ; )a 0 (2 
where 


transform variable of Laplace transform 
taken with respect to s Vt b 
a *= Laplace transformed of / and a 
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C(p) = transformed Theodorsen function 
k = wh V (reduced frequency 

ka = w,b/V 

7 = m mpb* (density ratio 

rae Xe —~ (Q/p 


and the remaining symbols are as conventionally 
described in the literature. 

In all the developments that follow the density ratio, 
uw will be assumed high enough (say uw > 15) to permit 
neglecting the virtual mass contribution to the A+ and 


aa terms, that is, 
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One manner of evaluating the flutter frequency and 
speed is to first determine the characteristic equation 
of Eqs. (1) and (2), next harmonic motion 
(p = tk), then determine the values of k and k, which 
make both the real and imaginary parts of the resulting 


assume 


equation equal to zero, 
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Low k Case 


For the conventional quasi-steady case, one assumes 


C(k = 
k — 0. 


first the imaginary, then the real part equal to zero 
results in the following expressions for the flutter 


frequency and flutter speed: 


a Yar + (wn We)*la — (1/2) ]Ja 
ee 2 
we V2 Fe + +a%—-U DE +e 
| 2) Q? yr - 4 4 © 
2 [r.? — O7(r,* — Fa°)] — (wr/we)*ra(1 — 2 








— 
Dita V bb \ 2,2? 


For the important case where w, w, is small (w, wy > 0 


, | ° 
J | ro 


burg V eB 


For the purpose of discussion, the first radical of 
Eq. (5) is designated A,, and the second radical, A». 
It is to be noted that A, corresponds, for high uy, 
identically with an empirical formula suggested by 
1940. Theodorsen and 


Garrick’s empirical formula was shown by comparison 


Theodorsen and Garrick? in 
with calculations to be a reasonable approximation to 
the flutter speed including the C(k) effect for the case of 
small w, wa, for the range of x, + a from —0.4 to 0.1, 
and for u > 10. The second radical, A», of Eq. (5) can 
therefore be looked upon as indicative of the effect of 
using quasi-steady rather than nonsteady air forces 
for the case where w, w, > 0. Fig. 1 gives a plot of 
the radical A, versus the significant parameter, t, + a, 
(roughly the c.g. location) for a typical value of rz = 
0.50. 


In Fig. 2, 


Conversion to other values of r, follows readily. 
the solid lines represent A, plotted versus 
Yr” = £7 


tz, + a for different values of 7, = Ss = 


a 


(roughly the nondimensional radius of gyration about 
the c.g.). 

From Fig. 2, it appears that the use of quasi-steady 
air forces results in lower flutter speeds than when using 
nonsteady air forces, the discrepancy getting worse as 
the c.g. approaches the mid-chord and as r, increases 


from When *, + a = O, the quasi-steady 


flutter speed indeed drops to zero. 


ZeCTO. 
For the range of 
0< % +a < 0.5, it can be shown by examination of 


the coefficients of the characteristic equation of Eqs. (1) 
and (2) and its discriminant that the system is unstable 


Voi 2) +a + Ita — (1 u) | V2 _ a + (z.. — a)- — (] ZA ae -+- g 


— (] u) | — (ay (es) 


a 


), the flutter speed reduces to 


(f, + a)* — (1/2)(% + a 


a 


For <=, + a> 0.5 


For values 


for all air speeds greater than zero. 
the flutter speed increases again from zero. 


of % + a below —0.5 + (1 yw), the flutter speed is 
infinite. 
The reduced frequency, k, associated with thes¢ 


quasi-steady flutter speeds is readily obtained from 
Eqs. (5) and (3). This varies from k = Oat % + a = 


at t, +-a=0. Itremains 


—0.5 + (1 uu) to k—> mt 
infinite for 0 < % + a< 0.5 and then decreases again 
as the flutter speed increases for %, + a > 0.5. It is 
interesting to note the great effect of the C(k) function 
on flutter of airfoils whose c.g. location is near the mid- 
chord but the minor effect when the c.g. is near the 
quarter chord. Another interesting phenomenon, that 
of flutter at zero air speed and infinite k when the c.g. 
the the 


chord, is also to be noted. The introduction of some 


lies between mid-chord and three-quarter 
additional damping to the system would raise this zero 


flutter speed to some small finite value. 


High k Case 


The possibility of low-speed, high & flutter evidenced 
in the previous case leads naturally to another ‘‘quasi- 
namely, setting C(k 0.5. 


steady’ approximation 


This is the asymptotic value for k > © and corresponds 
to the low-speed flutter examined previously by Biot 
and Arnold.’ Expansion of the determinant as before 


leads to 


Ww oe r (wy Ya la - (1 2)] 
- 2 = 0) 
i, Ve ae t+ le ta - 0) 
‘ ] = = > : * 
07[r.2 — 22(r,2 — £,2)] — (wa/we)ro2(1 — 2 - 
/ a 9 . = é ‘ 
bwa Vu Ver 2) +a+% — (1 w)] — (wn @,)*la + 1 2)] 
For small a, ‘wy, the flutter speed becomes 
: | Pe ; 
J | t ] 2(% +a — (1/2)}? 
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_ | | (S 
Vel > 404+K— Ww) We ee + le pa DP 


| which corresponds to the assumption that 
Expansion of the determinant and setting 
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The first radical in Eq. (S) is seen to be A, while the 
sues second radical is now designated A;. As in the case of 
n that ko, the quantity A; is a function only of t, + a andr 


and is shown plotted with the dashed lines in Fig. 2. 


It can be seen that the zero air speed, infinite & flutter 


setting 


O ZCTO 
itter now only occurs € 1e c.g. lies a 1e three-quarte depends rather on location of the c.g. 
flutt 7 h rs when tl g. 1 t the three-quarter lepends rather on location of tl 
chord point. As the region near this point corresponds 
to high k where the approximation that C(k) = 0.5 is Other Quasi-Steady Cases 
2 valid, the quasi-steady analysis gives an accurate To account for actual lift-curve slopes, the factor 
picture of the nonsteady flutter speed in this region C(k) is sometimes replaced by some constant, F, 
generally smaller than unity. For these cases, the 
flutter frequencies and speeds become 
w ran t+ (wn We)? [a 1/2)J,a — (1/2)[(1/F 1], 
=~ 9 = q 
Wa . Px —~F2+ (& +a)? — (1 2F)(% +a) + (1/4)[(1 F 
0? [ra” — 2(r.? — Fq°)] — (wn/wae)*ra(1 — 2 
F = 0 
j 0OWa VV M Verio 1/2) ta+Z— (1/p |— (w w,)*la + (1/2 ; 
and for small a, wg, 
| Re (t, + a)? — (1/2F)(% + 4 1 4)[(1 # 
> O.5 a | - : 1] 
bia V be 2/1 2)+a+2%,-—(1 w)] Fira? — ®" + (® +4)? — (1 2F)(& + a 1 4){(1 F 
alues 
ed is i o- —" 
The values of the modifying constant in Eq. (11 2 
: a For certain parameters, this divergence speed may 
—_ corresponding to A, and A; have not been plotted on , , 
1ese 6 ; back Nee ‘ very well fall below the flutter speed of the airfoil 
fees Fig. 2 for different values of Ff. However, it can be . 
aa shown that they will lie somewhere in between the a 
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the empirical Theodorsen-Garrick formula as ex 
pressed by A, was not intended to apply for t, + a 
It is interesting to note that 


the low-speed flutter of Biot 


much greater than 0.1 
for this case of low a, wg, 
and Arnold is independent of elastic axis position but 
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An Investigation of ‘Two-Dimensional 
Supersonic Base Pressures’ 


A. F. CHARWAT* anno J. K. YAKURA** 
University of California at Los Angeles 


SUMMARY 


An investigation of the base pressure behind wedges at Mach 
Numbers 2 and 3 in the laminar and the transitional regime is 


reported. Temperature and velocity traverses through the 


mixing zone are shown and exploratory investigations of the 
wake vortex by use of hot wires and flow-visualization techniques 
are described. It is found that the laminar two-dimensional base 


pressure agrees well with Chapman’s theoretical predictions 


The shear laver exhibits gross velocity distributions character 
istic of the free jet mixing zone, but also shows disturbances that 
originate in the expansion-turning of the oncoming boundary 
trailing vortex is observed, which is 


laver \n interesting 


explained in terms of nonuniform mixing rate in the wake 


INTRODUCTION 


— STRUCTURE of a wake at the base of a blunt body 
in supersonic flow is a particular case of a group of 
problems in which the field is determined by the inter 
action of a hyperbolic region with a region of dissipative 
flow. Another example of this class of problems is the 


she ck 


The particular importance of 


separated boundary layer encountered in 
boundary interactions. 
reaching a clear understanding not only of the phe- 
nomenological picture but also of its quantitative 
details stems from the characteristic sensitivity of the 
external supersonic flow to small changes in the shape 
of its boundary. The 


different flow fields makes this case an exceptionally 


fine balance between two so 
interesting one. 

From the practical viewpoint, the pressure in the 
wakes of blunt bodies is of considerable importance 
because it can account for as much as one third of the 
drag of medium-range supersonic vehicles. The corre 
sponding problem in two dimensions gained interest 
when it was shown that blunt trailing-edge airfoils can 
offer certain advantages. Besides the concern with the 
the flow in the subsonic wake 
the heat 


transfer rates and the cooling of jet nozzles located at 


base pressure proper, 


determines such critical phenomena as 


the base of missiles, and many associated problems. 
Accordingly, the supersonic wake (see Fig. 1) has 


received considerable attention. Early experiments 


led to an accumulation of data (summarized in refer 
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ences 1 and 2) which was found rather difficult to unify 
The first analysis which explained (at least qualita 
tively) this complex phenomenon was due to Croce 
and Lees.* Their formulation of the base flow structur 
emphasized for the first time the importance of th 
mixing process in determining the equilibrium between | 
the the Another 
approach was suggested by Chapman? and Korst’ and 


external flow and enclosed wake. 
also by Kirk,® but it is less general, although consider 
ably simpler. 

Crocco and Lees study the “‘free jet’ at the edge oj 
the supersonic flow downstream of the base. They use 
an integral approach whereby the essential momentum 
transfer is described by properly defined averages 
The ‘“‘dead air’ region behind the base is thus energized 
until it can support the recompression necessary to turn 
the supersonic stream back to the direction of the mean 
flow. No account is taken of the existence of fluid in 
the wake. 

The Chapman-Korst theory requires the conserva- 
tion of mass within the wake, where momentum transfer 
from the external stream induces a circulatory motion. 
In their the 
dynamics of the streamline which divides the two flows 


formulation these authors emphasize 
and express the equilibrium condition by stating that 
the total pressure along that streamline must just 
to the internal 


fiuid is returned and conserved within the wake while 


suffice support the recompression: 
the fluid external to that streamline escapes down- 
stream}. 

The aforementioned model leads to a simple expres- 
sion for the base pressure in two-dimensional flow, 


P» j l T iy | 21L1/ | ' 
P, i+ [(y — 1) 2] — U2) )S 


where P, is the pressure downstream of the recompres- 
sion region (equal, approximately, to the free-stream 
and Ll’, U/U, is the 


velocity at the dividing streamline just ahead of the 


pressure for slender bodies 


recompression referred to the external velocity. This, 
of course, does not constitute a solution because Ul’, is 


not known in general. However, some _ idealized 
models are available; Chapman uses his solution for 
laminar mixing supersonic jet free of boundary layer 
with stagnant air,’ a similar solution which yields 
U’, = 0.587, and Korst attacks the problem of turbulent 
mixing by use of an eddy-viscosity concept. Korst’s 
analysis is formulated in general, but carried out only 


for no oncoming boundary layer. 
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Schlieren photographs of the base wake. (A 
R, = 610,000; (D) MW 


Both analyses assume an isoenergetic and isobaric 
mixing process. 

There seems little doubt that the outline of the flow 
field behind blunt-body bases described above is con- 
However, there are many details 
which To with, the two- 


dimensional laminar wake which is solved most com- 


ceptually correct. 


remain obscure. begin 


pletely within the framework of these simplified the- 
ories has never been subject to experiment;* only the 
axisymmetric flow 


corresponding three-dimensional 


was studied in the laboratory.*:* Furthermore, very 
little attention was devoted to discovering how the 
actual velocity distributions across the mixing layer 
The only data 
10, 


compare with the simplified models. 
to the 


which became 


known authors was obtained in reference 
available during the course of the present 
investigations. This information is necessary to guide 


further work and refine the understanding of such 
problems as the effect of boundary layers on the struc 
ture of The experi 


ments reported herein were designed to investigate some 


the inclosed ‘‘dead air’’ region. 


of these details. 


BASE PRESSURE TESTS ON WEDGES 


The base pressure measurements were conducted in a 
1.98 and 17) The 


models were geometrically similar wedges with a 4:1 


3 by 3 in. tunnel at 1/7, 2.78. 


chord to base height ratio. The variation in Reynolds 
Number was obtained partly by varying the chord 
length of the models and partly by throttling the tun 


nel. The largest model had a 2-in. chord and an aspect 


while the smallest chord was 0.3125 in. 
All models were canti 


ratio of 1.5, 
with an aspect ratio of 9.6. 
levered from one of the tunnel walls and spanned its 


* Reference 4 reports data for a backward facing step. How 
ever, there exists some evidence that the base pressure is higher 
in that case than in the corresponding free wake 
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M=28,R 








155,000: (C) AT = 2, 


= 580,000 B) VW =2.8. R 


= 2 R,. = 130,000 


STREAKLINES ON GLASS SIDE WALL 
ii INITIAL LAMPBLACK SMEAR 





IN M* 2 TUNNEL 


2 INCH WEDGE 
Fic. lb. The model installed. Tunnel side wall raised show 
ing turning of the flow and the recompression shock by lamp 
black streaks. Note interaction of ompression shock with 
tunnel floor boundary laver 


rec 


width. The pressure taps in the model were brought 


out through the side wall at the model support. The 
entire assembly was two-dimensional. 
Preliminary tests were conducted to verify the 


uniformity of the flow and the pressure along the span 
of the two largest models. In both cases negligible 
differences existed along the entire length except very 
close to the wall. Both the base pressure and the sur 
face pressure on the wedge just ahead of the base were 
checked. Additional tests 


that the boundary layer on the surface of the models 


flow-visualization showed 


was laminar and uniform throughout the range of in 
vestigation except at the tunnel wall; however, some 
interesting three-dimensional disturbar:ces in the wake 
These are discussed in the last para 


the 


were observed. 


Further measurements of base pressure 


along the height of the bases were made, showing, in all 


graph. 








124 JOURNAL OF THE 
cases, that it is uniform both in the turbulent and in the 
laminar regime. 

The measurements of base pressure are reported in 
Fig.2. The present data is extended into the turbulent 
regime by measurements of Gadd et al.'! who used, at 
M = 3, M=2 
Gadd’s data is for a more slender wedge (ch = 8). 
Noting the fairly well substantiated correlation of base 
pressure in the transitional regime with the parameter 
5h = chvVR, (see reference 4), it would appear that 
Gadd's data is displaced somewhat towards lower 


the identical wedge geometry. At 


Reynolds Numbers in comparison with ours, such as 
would occur if the turbulence level were larger in their 
tunnel (see Fig. +). Nonetheless, the agreement is 
good. For comparison, curves for a cone-cylinder*® are 
also shown. 

Without exception, all the evidence indicates the 
familiar rise in base pressure ratio from a (constant) 
“turbulent wake” value to a constant higher plateau for 
Reynolds Numbers below critical. This rise is ex- 
plained by the motion of the transition to turbulence in 
the free jet, and the consequent variation in momentum 
transfer to the shear layer; below the critical Reynolds 
Number the transition point is presumably located 
downstream of the recompression region. This is the 
region to which Chapman's theory should apply most 
closely. 

The scatter in the present measurements is small (at 
least two distinct tests are reported for each model) but 
there is a recurrent difference between the _ base- 
pressure ratios on different models at a fixed Reynolds 
Number. This difference is due, in part, to a slight 
change in Mach Number with Reynolds Number. It 
varied from 1.99 to 1.96 and from 2.78 to 2.66, respec- 
tively. A correction [using Eq. (1)] was tried, but it is 
too small to compensate for the observed differences. 
Special precautions had been taken to prevent non- 
uniformities in surface finish or leading-edge bluntness 
among the models. It seems justifiable to explain the 
discrepancy by a boundary-layer effect, the thickness of 
which (at a fixed length-Reynolds Number) is not 
reduced similarly to the model geometry. 
the appropriate boundary-layer 
The thickness proper 


The question of 
correction remains uncertain. 
would appear to describe the oncoming velocity gradient 
and its effect on the expansion-turning of the flow at the 
edge of the base; to describe the increased resistance to 
momentum transfer from the mean flow, the momen- 
tum thickness would be appropriate. The 
smoothing of disturbances in the profile along the mix- 
ing layer is function of Y and +~/X for laminar and 
turbulent jets respectively, where Y is the distance 
The base height (/:) can be taken to be 
In the transition regime 


more 


along the wake. 
roughly proportional to LX. 
Chapman's parameter 6/4, which correlates the data 
well, can be interpreted to represent a stability criterion 
for the boundary layer (function of its thickness) related 
to the distance along the shear layer available for 
transition. It does not follow that the same parameter 
should apply to the fully developed turbulent regime or 
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the pure laminar one. Fuller et al.'’ use a_ linear 


extrapolation of the momentum thickness to bas 
height ratio to correct their data (turbulent wake) and 
compare it with Chapman’s;'! Gadd,'! who used much 
smaller models for which the momentum thickness is 
essentially zero without correction, checks the afore- 
mentioned results very closely. It appears that for 
turbulent wakes the ratio 6 is appropriate. The 
present data, in the laminar regime, is shown related ti 
the parameters in question on Fig. 3. It is evident that 
the best correlation is obtained against the boundary- 
thickness Moreover, the theoretical 
value of U, = 
the experimental curve. 

A set of base-pressure measurements was obtained 


proper. 


layer 
0.587 for laminar mixing layer falls on 


with the leading edge of the models roughened by 
attaching to it a strip of 00 emery cloth. 

The roughened wedge surface tests were conducted 
They are 


on the smallest and the largest models. 


shown in Fig. 4. The effect on the base-pressure vari- 
ation was to shift the curve towards lower Reynolds 
Numbers by increasing the tendency for turbulent 
transition. The became 
higher in accordance with the aforementioned bound- 
The turbulent wake flow regime 


seemed to be less affected. 


final laminar base pressure 


ary-layer effect. 


THE MIXING LAYER 


The key to the base-wake problem lies in the dy- 
namics of the free mixing layer which separates the 
external flow from the internal vortex. It is necessary 
to know its mass and momentum distribution and its 
stability characteristics, as well as the influence of the 
oncoming boundary-layer parameters on its develop- 
ment. 

Fig. 5 shows a number of successive traverses of the 
shear layer at two Mach Numbers and two Reynolds 
Numbers. The were made with a 
flattened total-head tube (overall thickness 0.004 in. 


measurements 
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which was lined up with the direction of the local flow 


in the shear layer by visual comparison on the schlieren 
screen. The data, which are reported in terms of local 
Mach Numbers, were computed on the assumption 
that the static pressure was constant and equal to the 
base pressure. This is meaningful only up to the lip 
shock or the base expansion fan. The structure of the 
hyperbolic flow, as indicated in Fig. 5, is selfconsistent, 
and the location of the measured dip in Mach Number 
at the outer edge of the layer coincides with the lip 
shock as viewed on the schlieren. 

It is immediately evident that in the early stages of 
development, the shear layer resembles the theoretical 
profile only in its gross characteristics. In particular, 
a ‘cusp’ in the velocity distribution is apparent. It 
diffuses but persists even across the recompression zone. 
This discontinuity is much stronger at the higher Mach 
Number. It is interpreted as follows. The oncoming 
supersonic flow is rotational and, in undergoing the 
expansion-turning at the base, a system of waves of 
both families is generated. When the Mach Number is 
greater than 7/2, that is, in the outer portions of the 
boundary layer, an expansion is reflected as an expan- 
sion into the slower layer and a compression wave is 
emitted into the faster one (away from the base). 
These shocklets coalesce to form the well known lip 
shock. However, when the Mach Number is less than 
V 2, that is, in the inner portion of the rotational layer, 
an expansion wave is reflected inwards as a shock. The 
net effect of the waves reflected off the 17 = +/2 layer 
(not a streamline) during the expansion is to form the 
slip-line noticed in the measurements. 
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Fig. 6 shows some profiles before recompression 
similarity parameter 
0.65 at y = 0, 


plotted against the turbulent 


The origin was chosen so as to make Ll’ 


(For clarity, a limited number of specific points was 


plotted instead of the entire profile.) The value of the 


scale factor o was chosen so as to make the inner 


boundary of the mixing region coincide approximately 
with that derived theoretically. This value is o 
20.0 for all three profiles shown in Fig. 5. 

The incompressible mixing layer was found to follow 
with o between 12 and 14 


the theoretical solution 


Investigations of mixing!*® of compressible fluids at 


M = 


between the incompressible and compressible profiles 


1.6 proved that there is no essential difference 


particularly in the inner low-speed portion of the layer. 
The scale factor o was found to be higher, approxi 
The present magnitude of ¢ = 20 does not 
Note that 
the aforemen- 


mately 15. 
seem to constitute a serious difference. 
according to the mixing layer theory 
tioned factor is related to the ‘‘mixing length” as the 
inverse 3 2 power. Consequently a higher value of ¢ 
indicates a lower turbulence level; this may be associ 
ated with the increased stability of supersonic jets, or 
more directly, with the fact that the Reynolds Number 
at which these profiles were measured was not yet high 
enough to place the flow in the “fully turbulent 
regime. 

For comparison, the profiles close to the base were 
also investigated. They were found to depart sub- 
stantially from both the laminar and the turbulent 
seems 


distribution. A graphical comparison 


fluous; it suffices to state that the profile has a thickness 


super- 


which substantially exceeds the turbulent and therefore 
also the laminar value. To fit it to theory in terms of 
the turbulent similarity parameter one would choose a 
value of o in the neighborhood of 4. This fact can also 
At AJ 
escape interpretation because of the proximity of the 


expansion fan, but they are well defined at JJ = 2; 


be seen by inspection of Fig. 4. 3 the profiles 


the shear layer does not grow as rapidly as would be 
expected. 

Finally, it is worthy of note that the inflection in the 
boundary-layer profiles denoted before as the internal 
slipline finds itself located quite near the value of lx = 
U/U, 


Since the dynamic energy of that streamline determines 


which corresponds to the ‘“‘dividing-streamline.” 


the recompression, the occurrence of this disturbance 

may bear on the resulting base pressure. 
Unfortunately, we were unable to study the laminar 

mixing region profiles because of the small size of the 


model. 


FLOW IN THE WAKE 


The mass induced by the mixing is small, and there- 
fore the velocities and energy levels in the recirculating 
flow are low and do not influence the dynamics of the 
field. The 
separate problem determined entirely by the conditions 


subsonic ‘‘base vortex’’ constitutes a 


at the boundary. An elegant analysis of flows with 
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closed streamlines was published recently and 
Batchelor'® proposed a model of laminar wakes based on 
this theory. It predicts that the vorticity is constant 


in the region, that is, if a single stationary vortex 
existed there, it would be of the solid body type. 

We conducted exploratory hot-wire traverses of the 
wake region which are given in Fig. 7 as ratios of the 
heat-transfer coefficient \ to its value at the centerline. 
We did not feel justified in reducing the data to veloci- 
ties because of the uncertainty in the wire calibration. 
However, the relation between \ and Ll’ is monotonic, 
ind the increase in A can be interpreted as an increase in 

Thus no stagnation surface dividing the down 
stream flow from the reverse flow in the wake 1s evident 
from the measurement. The reason is that the meas 


urements represent an average root-mean-squared 


elocity because the wake flow is unsteady. Observa- 
tions near the inner edge of the mixing layer showed 
large aperiodic fluctuations even at the lower Reynolds 
Number. These waves of the boundary are imposed 
m the wake; their period is long with respect to the 
that is, they do 


turbulence in the 


characteristic time of the shear layer 


not necessarily represent mixing 
layer) but short with respect to the slower wake flow 
ind they cause it to break up. An unstable vortex can 
observations of 


the 


iso be inferred from an aperiodic 


vortex shedding into wake observed on spark- 


schlieren photographs by Chapman.' This situation 
probably persists to very low Reynolds Numbers, much 
below the critical, where the shear-layer proper is 
laminar. 

It seemed interesting to verily the existence of the 
recirculation’ of mass in the wake. The base of the 
wedge was covered with a strip of a pH indicator and 
ammonia was admitted into the wake through a total 
When the tube mouth was brought 
“throat,” the 
carried the vapors to the base as evidenced by a change 
the The effect was roughly 


uniform over the base height; there was no time delay 


pressure tube. 


forward past the wake recirculation 


in color of indicator. 
in color change near the edge as compared to the center- 
line. This also points to a high degree of turbulent 
Incidentally, the region 
near the ‘‘throat”’ the 
would not penetrate into the wake at all coincided with 


mass-mixing in the wake. 


downstream of which vapors 
the recompression region as viewed on the schlieren; 
thus the wake was closed. 

Fig. 7 also shows a total temperature traverse of the 
The wake 
This is high enough to justify the assumption of 
isoenergetic flow. The fact that the recovery tempera 
ture in the base wake is high and, also that the convec 


wake. factor was found to be 


0.945. 


recovery 


tive heat-transfer characteristics are poor suggests that 
the base region may lead to aerodynamic heating 
problems at hypersonic speeds which are at least as 
severe as in the forward stagnation region. 

Particular attention was given to the study of the 
uniformity of the pressure along the height of the base. 
No significant pressure distribution was found in either 


the turbulent or the laminar regime. The same was 
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true for the pressure along the centerline in the stream 


course, near the recompression 


direction 
region (incidentally, the data is not accurate there since 


except, ol 


significant disturbances are introduced by the probe 


The two-dimensional wake is thus isobaric, which is not 


so in the case of bodies of revolution It would appear 


that the gradual recompression characteristic of con 
verging conical flows, which is responsible for the higher 
rise to the 


mean value cf base also gives 


radial distribution of pressure in axisymmetric wakes, 


pressure, 


OBSERVATIONS OF NONUNIFORM WAKE MIXING 


An interesting observation was made during at 


tempts to visualize the flow pattern. Fig. S shows the 
streaklines over the base; these are lampblac k deposits 
obtained by coating the base with a kerosene-lampblack 
axis parallel 


suspension. <A trailing vortex with the 


to the flow is seen to be shed from the wake For a 
short time before the kerosene evaporated we could 
Note that there are 


A pair can be faintly 


actually observe the circulation 
actually 4 pairs of these vortices. 
rhe other two pairs, of 


seen at either end of the span. 


BASE 


LEADING 
EDGE 





Lampblack streaks evidencing the shedding of a vortex 


downstream from the base; JJ = 2, R, = 600,000 
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which only one was properly illuminated and registered 
clearly on the photograph, are located symmetrically on 
each side of mid-span. 

The phenomenon, as well as the location of the 
vortices, was consistently reproducible at J/) = 2. 
However, at higher Mach Number we were unable to 
observe any consistent flow structure whatever over the 
base. (It is possible that the violence of the starting 
shock and the higher shear rates tending to carry the 
fluid downstream prevented a successful application of 
this technique. 

The trailing vortex is a three-dimensional disturbance 
with its axis perpendicular to the mixing vortex; 
incidentally, the latter was also observed by lampblack 
streaklines formed on the tunnel side walls. 

The shedding of vorticity oecurred at the intersection 
of the base with the Mach line originating at the junc- 
tion of the leading edge with the wall. It seemed that 
this slight disturbance (which did not affect the surface 
boundary layer or modify the surface pressure notice- 
ably) was sufficient to alter the stability characteristics 
of the free-shear layer and induce an earlier transition. 
Thus the mixing in the free-shear layer was increased 
and resulted in a stronger (mixing) vortex in the 
affected portion of the wake span. The excess vorticity 
had to be shed at the junction of the disturbed and the 
undisturbed regions to satisfy vorticity conservation 


laws. This gave rise to the appearance of the trailing 
vortex. To verify this hypothesis we glued sandpaper 


strips over arbitrary portions of the span to artificially 
destabilize the boundary layer. Indeed, we found that 
a vortex was shed downstream at the edge of the dis- 
turbed region every time. 

Despite the existence of this three-dimensional dis- 
turbance, the base pressure remained uniform across 
the span within the accuracy of the readings. The 
three-dimensional motion appeared to be too slow to 
affect the basic wake pattern. Still, to ascertain 
definitely that the pressures, and, principally, the 
critical Reynolds Number as determined by the stabil- 
itv of the mixing layer, were representative of truly 
two-dimensional flows, two thin plates were brought in 
from downstream against the base. These end plates 
(perpendicular to the span) effectively isolated the 
central undisturbed portion of the wake from the dis- 
turbed end-zones. The entire series of tests was re- 
peated and it was found that neither the magnitude of 
the base pressure nor the critical Reynolds Number 
were altered. 

It was concluded that the data presented in this paper 
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were not vitiated by the existence of the three-dimer 
sional trailing vortex. However, it is possible that 
such three-dimensional disturbances do affect transitio; 
in the wake over the local portion of the span wher 
they occur. In our tests, the model used in the critical 
Reynolds Number region had so large an aspect rati 
that the disturbance may not have diffused to the center 


of the span where the pressure tap was located. 
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A Remark on Stability Equations 


Erich Sponder 

Project Engineer, Bolkow-Entwicklur 
Germany 

September 13, 1957 


ger tuttgeart rport 


US CONSIDER the simpiified case of a purely vawing air 


| EI 
plane, as shown in Fig. 1. The 


constant without 


tirplane may proceed hori 


zontally at speed experiencing any drag, 


thrust, or side force but will experience a yawing moment .N if 


there is generated an angle of sideslip 8 by a small disturbance 


\ damping in yaw, on the contrary, will not exist. Consequently, 


the flight path of the c.g. must be a straight line along which 


vawing oscillations may occur since the airplane is supposed to 


be statically stable: ON/Os 0 (weathercock stability With 
a moment of inertia / about the ¢.g., the equation for the yvawing 
motion becomes 
I(d*3 /dt 0.V/03)3 i 
which describes an undamped permanent oscillation 
Now, let us add a side force, or cross-wind force, J}, as indi 
cated in Fig. 2. After a disturbance in vaw, the side force } 


0}'/03)8 (where 0)'/08 > O) pushes the airplane off its course 
to the right, so that there will be a small departure from the orig 
The the 


ensuing motion now has to distinguish carefully 


inal direction of the e.g. (dashed line description of 


umong the three 


angles 8, x, and ¥, where ¥ = 8 + x \gain, a damping in vaw 


will be left out Then we have 


I(d*y dt?) — (0.V/03)3 i) 


different from the former expression, and, with the mass m of 


the airplane and its speed JV, 


mV (dx /dt 0}'/03)3 


of the first equation as /{(d*3 


can eliminate d?x ‘df? by 


Writing the term /(d*y /df? 


dt?) + d*x/dt*|, we differentiating the 





Contributors Please Note 


Readers’ Forum items must be confined to the 
equivalent of one page in the Journal. Contribu- 
tions that exceed this limit will be returned to the 
authors for condensation and rewrite. To avoid incon- 
venience and delay, please adhere to the following 
specifications: 

(1) Five, double-spaced, typewritten, manu- 

script pages (8'/, by 11 in.), including 
wide margins, formulas, and headings 
equa! one printed page. 
For every illustration, deduct at least 
one-half or as much as one manuscript 
page, depending on the size of the 
illustration. 


(2 


~— 


























JOURNAL OF THE AE 


RONAUTICAL 





Flight Path 





Fic. 1. A purely yvawing airplane along a rigid, straight flight 


path 








Original Direction 


Fic. 2 vawing airplane along a yielding, curved 


flight path 


A purely 


second equation; hence 


I(d*8/dt?) + (OY/ogs) (I/mV) (dB/dt) — (ON/OB)B = O 


This result is quite remarkable. Although no damping in yaw 
was introduced at all, the 6 oscillation of the airplane about its 
curved flight path will gradually disappear since the coefficient 
of d8/dt is essentially positive 

This is the formal result of a treatment which proceeds from 
simplified to more complicated modes of motions as is usual in 
lateral dynamics. Yet, it doesn’t help to give a clear physical 
understanding of this particular damping effect. In order to 
shed some light on this question, one may try to compare the 
curved flight path of the yvawing airplane with the motion of a 
without friction along a 


solid body which slides horizontally 


curved wall. But such attempts cannot yield any similar re 


sults, because no energy is dissipated 
The correct explanation of this damping effect seems to be 


quite obvious when it is found by closer examination. The side 
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force ¥ can be produced only if there is an equivalent amour 





momentum imparted to the surrounding air in the opposit 


direction 


away some energy which clears up the decrease of the yawing 


oscillation 

\ complete treatment of the airplane dynamics contaii 
course, terms of the same significance, though in a more involve 
\n inspection of the stability quartics shows them dis 
Aerodynamic T) 


IS6 for t 


form 


tinctly in the coefficient B (see Durand, 
vol. V, p 


lateral motion 


171 for the longitudinal motion, and p 


The shortened notations s, and y, mean exact 


the slopes of the lift curve and of the side-force curve respectivel 


Besides, Fig. 1 corresponds to an experimental investigation 


a wind-tunnel 


the vawing mode by test, while Fig. 2, wit 


more lateral degree of freedom, better approaches the airplane's 
motion in free flight Since this effect is of the same nature 
in induced drag, which is necessarily connected with a lift 


sideslipping, pitching, or reeling motion of airplanes, projectiles 


ind the like, will experience a damping which sometimes may by 


considered as a desirable gift for some improvement of the d 


namic stability ‘ 


An Approximation of the Boundary of a 
Supersonic Axisymmetric Jet Exhausting 
Into a Supersonic Stream 

E.S. Love 

Head, Hypersonic Tunnels Branch. Langley Aeror 

NACA, Langley Field, Va 

September 20, 1957 


Laboratory 


SYMBOLS 
V/ free-stream Mach Number 
Vv jet Mach Number immediately ahead of exit 


3 boat-tail angle 


semidivergence angle of conical nozzle 


radius of jet at exit 

distance normal to jet axis 

distance from plane of jet exit measured parallel to jet axis 
jet static pressure immediately ahead of jet exit 
free-stream static pressure 


ratio of specific heat 


INTRODUCTION 


OR SHORT DISTANCES DOWNSTREAM Of the exit of an axisym 

metric supersonic jet exhausting into a supersonic stream 
i.e., distances of the order of one jet exit radius—the boundary) 
between the jet and the stream may be satisfactorily predictec 
by the use of the boundary of the jet that exhausts into still air 
at a pressure ratio which produces the same initial inclination 
of the boundary as the jet exhausting into the supersonic stream 
However, for appreciable distances downstream of the jet exit, a 
satisfactory prediction of the boundary necessitates, at present 
a tedious application of the method of characteristics. The need 
for a simpler solution stimulated the examination of an approx! 


mate method summarized in this note 


DISCUSSION 


The approximate method is proposed for flows in which the 
local Mach 
ciently high to permit use of shock-expansion techniques. The 


Numbers on either side of the boundary are suffi 
method is restricted to the portion of the boundary ahead of tht 
maximum diameter of the free jet (the method falls apart at th 


maximum diameter and beyond), to positive initial inclinations 


of the boundary, and, probably, to moderate static pressure 
ratios. Both the expansive turning and the compressive turning 
of flow adjacent to the boundary that occurs after the boundary 


has reached its initial inclination are assumed to be isentropi 


The problem of calculating the boundary shape is simplified by 


And this ‘‘sidewash,”’ as it is named in Fig. 2, carries 





| 








Fic 





READERS 


sidering the shape to be determined by only three effects. 


umount These al an expansion of the ambient flow by turning along 
opposite F ye boundary, a compression of the jet flow by turning along the 
~, Carries undary, and an effective expansion of the overall jet flow in 
le Vawing | poceeding downstream that is governed by the change in cross 
tion irea of the jet flow The latter two effects are coun 
ntains, of | racting but produce a resulting effect which, within this ap- 
involy simation, must equal the expansive effect of the ambient 
them d 
tC Ti The initial calculations follow a conventional characteristics 
OD Tor t utiol the initial inclination of the boundary is calculated 
n exact nd assumed to hold over a short distance downstream of the jet 
ectivel <it (for 0.1 r; in these calculations), and the initial properties of 
§ation of F she flow on either side of the boundary are determined \t this 
wit int the approximate method is introduced. It is a step-by 


ir plane's tep procedure in which a small change in the slope of the bound 


lature as rv is introduced and the values of x/r; and satisfving this 
lift, ar inge in slope by means of the aforementioned effects are cal 
ojectiles ulated from one-dimensional relations and Prandtl-Mever 
> may be F turning \ standard change-in-slope increment of 1/2” was used 
the d the examples to follow; other calculations showed that 1 


nerements changed the approximate solution negligibly, while 
nerements were unsatisfactory Details of the calculations 
vill be omitted; briefly, a typical step in the calculations involves 
calculation of the change in pressure of the ambient flow at 


1 expanding 1/2°, the change in pressure of tl 


boundary 1 


e 
t flow at the boundary in compressing 1/2°, and the change 1n 
rea of the jet flow (expansion of the jet flow) that is necessary 

1g to create a pressure balance at the boundary From this newly 

etermined area and the slope of the boundary, the values of 
ind may be calculated Throughout the calculations 
effective throat area that is used to determine the values of 

s that given by the Mach Number of the jet flow at the point 

here the approximate method is first applied (at a = ().1 in 

ese calculations The use of conventional flow tables simpli 

s the calculations 

Comparisons of the theoretical exact boundary, given by the 

hod of characteristics, and the approximate boundary art 
wn in Fig. 1 The exact boundary for Fig. 1(b) was obtained 
reference 1; the value of 7 for the jet flow and ambient flow 

f this calculation was 1.405. For the remaining boundaries, 

s 1.400 rhis small difference in 


For Fig. l(a) the Mach Numbers of the jet flow 


is of no consequence in the 
is mparisons 

inbient flow adjacent to the boundary just downstream of 
For Fig. 1(b) these Mach 


In view of the lower 


e exit are 3.30 and 3.57, re spectiv ely 
Numbers are 2.62 and 1.72, respectively 
Mach Numbers involved, it is not surprising to find that the 


greement between the exact and the approximate boundary is 


seis fe oorer for Fig. 1(b). To facilitate this comparison, portions of 
ream 


undary 
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: ; Fic. 1. Comparison of exact boundary with approximate 
ropic boundary. 


ed by 
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the boundaries of Fig. l(a) have been presented in Fig. I(c) t 
the same scale as those in Fig. 1(b). 

In spite of the gross simplifications involved, the approximate 
method yields a fair prediction that may be useful in engineering 
applications. A manual calculation of the approximate boundary 
requires about 1 hour (with the aid of flow tables) as compared 
to several days for a manual calculation of the exact boundary 
In all cases examined, the approximate boundary falls beneath 


the exact boundary 
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Estimating Aerodynamic Characteristic Times 
in Hypersonic Flow 
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O~ OF THE MAJOR PROBLEMS associated with the flight of an 
f the 


object at hypersonic speed is the determination o 
effect ol the shock wave on the air Su) that the state of the dis ‘) 
ciated and ionized gas may be evaluated in the vicinity of the 
body 

If one wishes to determine whether thermochemical equilib 
rium exists at a particular region in the flow field, for example, 
at the outer edge of the boundary laver, not only must the rates 
be known at which an equilibrium state is approached, but the 
aerodynamic characteristic times must also be known for eacl 
point in the region Phe aerodynamic characteristic time is de 
fined as the time required for a gas particle to pass through the 


shock and reach a particular point in the flow field 


rhe presence of a blunt body completely distorts the velocity 
field downstream from the shock so that the aerodynamic char 
acteristic time is a function of the particle path selected and does 
not vary linearly with distance Hence, the treatment in refer 
ence 1 cannot be used directly to determine the departure from 
equilibrium at the outer edge of the boundary laver Note that 
the latter must be known if one is to apply the outer boundary 


conditions properly to a heat-transfer analysis 
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This note concerns itself with some first estimates of the aero 
dynamic characteristic times in the vicinity of the stagnation 
point. In Fig. 1 are shown nominal boundary-layer thicknesses, 
computed for a 2-ft. radius sphere. These thicknesses were com 
puted over a wide range of altitudes and Mach Numbers on the 
assumption that chemical reactions are frozen and that the flow 
is laminar.* It appears that the nominal boundary-layer thick- 
ness is a sensitive function of altitude but is relatively insensitive 
to Mach Number. 

For the sake of simplicity, the thermochemical problem is split 
from the gas-dynamics problem. It is further assumed that the 
structure of the velocity field, at least at the stagnation point, may 
be determined from a hypersonic approximation Note that the 
detachment distance has been assumed to be given by 6/Rg = 
0.13, which seems to be in agreement with recent theory and ex- 
periment.‘ Fig. 2 then follows directly. 

For the characteristic times indicated, it appears that there 
may be insufficient time for thermochemical equilibrium to be 


established at the stagnation point at the higher altitudes 
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On the Conservativeness of Various Distributed 
Force Systems 
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-_—- NECESSITY of determining the conservativeness of force 
systems in instability problems is discussed in reference 1 


AERONAUTICAL 


SCIENCES FEBRUARY, 1958 


It is shown that, whereas kinetic methods are generally app! 
for the determination of instability loads, the statical metho 
usually employed are valid only for conservative and nongyr 


scopic systems. Small changes in the character of the loading 


could make an otherwise conservative system nonconservatiy¢ 
and cause a large change in the magnitude of the buckling loa 
The buckling load of the cantilever column example in referenc 
1 is, for tangential end loading, eight times that for conservative 
constant directional loading 

Two papers have appeared recently* * which have treated 1 
buckling of circular rings (or infinitely long cylindrical shells 
under various distributed force systems. Different values of the 
instability load were obtained for the various systems. Since 
there appears to be only slight differences in the character of thes 
systems, it is of interest to examine them in regard to their ¢ 
servativeness and their properties with respect to the ring buck 
hydrostatic pres 


ling problem. These force systems are: l 


sure, (2) constant directional pressure, and (3) centrally directed 


pressure. The magnitude of the pressure is assumed to be con- 


stant during buckling. Since a circular ring buckles in an inex 
tensional mode, the magnitude of the total applied force also re 


mains constant. The direction of the resultant force on each ring 


» 


element does change for systems 1 and 3. It is this change i 


force direction that accounts for the difference in buckling loads 


between the systems 


HYDROSTATIC PRESSURE 


The work done by a constant hydrostatic pressure p acting 
upon an enclosed volume as the volume contracts by an amount 
AV is pAV 
such a system is conservative appears in texts on hydrodynamics 


} 


A proof of this expression as well as a proof that 


e.g., pp. 13, 14 of reference 4+. This system is also conserva 
tive for variable pressure if the pressure can be expressed as a con 
tinuous and single-valued function of the volume 

The buckling load of a circular ring under hydrostatic pressure 


is, in usual notation, 
g=3EI/R l 

where g is the load per unit circumferential length. The work 
done by the hydrostatic force system during buckling involves 
linear and higher terms in the buckling displacements and their 
derivatives. The linear terms must cancel out in the expression 
for the total energy since the buckled mode is an equilibrium 
position, and it is the quadratic terms that enter into the buckling 
criterion. The quadratic terms in the work done by the hydro 
static force system are given on p. 142 of reference 2 

On the basis of physical reasoning, it seems that the work done 
by the hydrostatic force system during ring buckling would be 
greater than that of any other conservative system. (A system 
in which the force is defined to be in the direction of the instan- 
taneous motion of the element on which it acts would depend 
upon the path and hence be nonconservative The ring buckling 
load for hydrostatic pressure would therefore be expected to be 
lower than that for other uniformly distributed conservative forcé 
systems 


CONSTANT DIRECTIONAL PRESSURE 


The constant directional pressure force system is defined as 
one in which the direction of the pressure remains constant ir- 
respective of any changes in the position of the elements. Since 
the work done by a force that is constant in direction and magni- 
tude depends only upon the initial and final configurations, this 
system is conservative. Although such a force field is easily de 
fined, it would be difficult to devise an experimental setup in 
which such a field is applied to a ring 

During ring buckling due to this force system the force vector on 
each ring element translates but does not rotate. The expression 
for the work done therefore involves only linear terms in the dis 
placements and their derivatives since the quadratic terms aris¢ 
These linear terms, as stated 


Buckling in 


due to rotation of the elements. 
previously, do not enter into the buckling criterion 


READER 


applicab| 


il method this ca is in the Euler column, is governed only by the energy 
1 nongyr f the induced bending stresses and the part of the work done by 
he loading the prebuckling membrane stresses which involves nonlinear dis 
MSETV ative nlacement terms 
kling loa The instability load of a circular ring subjected to a uniform 
1 refereng circumferential pressure of the kind considered here is given by 
Iservative Eq. (24) of reference 2 (upon setting vy = 0) as 

g = 4EI/R 2 
reated t 
cal snel 


CENTRALLY DIRECTED PRESSURI 
ues Of the 


ns. Sir This force system has been used in reference 3 in relation to the 


ice 
er of the ring buckling problem. During buckling the force vector on each 
their cor ring element remains directed toward a common pole (the ring 


The work done 


i unit force during buckling is the product of the force 


nter ind therefore both translates and rotates 


ing buck 


atic pres. by su 
> directed i the change in the distance of the element from the pole Phe 
o be con- work is therefore independent of the path so that the force field 


in inex. | is conser itive If the force intensity were not constant but 


e also re. lepended in a continuous and single valued manner upon the dis 
each ring placement from the pole, then the svstem would still be conserva 
ach ring 

-hange in tive—e.g., the general central force problem 

ing loads The buckling load for centrally directed pressure is given by 


> 


Eq 37) of reference 3 as 
g = +.5 EI/R 3 
The part of the work done by this force system which involves 
p acting quadratic displacement terms is negative for the ring buckling 
amount 


‘oof that 


lisplacements Eq. (15) of reference 3 is actually a potential 
nergy expression and hence is the negative of the corresponding 
[vnamics work term The buckling load for this case can therefore be 


onserva expected to be greater than that of case 2. It is of interest to 


iS a con 


note that the quadratic terms in the expression for the work donc 


by the external force system are positive for case 1, zero for case 


pressure ind negative for case 3 
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A New Method for Computing Drag 


syste . 6 . . 
= aa Coefficients from Ballistic Range Data 
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depend ee i 
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p« SENTLY available methods for computing the drag coefficient 


from time-distance data obtained in ballistics ranges are 





1ed as satisfactory but laborious. The drag foree cannot be regarded as 
int ir- constant over long intervals of model flight because of the vari 
Since ition in dynamic pressure with model velocity. The resulting 
nagni- deviation from uniformly decelerated motion has led to one pro 
s, this cedure in which the time-distance record is represented by a power 
lv de series expansion which is fitted to the data by least squares! and 
up in mother in which iterations are performed starting with uni 
formly decelerated motion as a first approximation.’ 
tor on Some of the author’s associates were using a graphical procedure 
eSsi0n lor a first check on time-distance data from the Ames supersonic 


e dis lree-flight wind tunnel \verage velocities in successive intervals 


arise | were plotted as a function of distance to determine the smooth 
tated ness of the data. Some of the curves were remarkably linear, 
ing in ind this suggested that one should simply take the slope of the 
] 


line as a measure of the drag coefficient However, the equation 
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of decelerating motion 
Cpo(1/2)pV?A = —n dV /dt = —mV(dV/dx ] 


shows that dV'/dx will not be constant except in the limit of 
vanishing velocity loss, so this procedure will have shortcomings 
for cases where the velocity loss is verv large It should be noted 
that the variations in dV /dx are about half as great as the vari 
ations in dV/dt so that the plot of velocity against distance does 
have merit compared to the plot of velocity against time 

These observations caused the author to look for some pres 
entation of the data that would be linear even in the case of 
large velocity loss and such a presentation was quickly found 


By expressing Eq. (1) as 


where K 


it can be seen that /m V is linear in x as long as Cp ts constant 
Thus, a plot of /n V against x should give a straight line whose 
slope is — A Cp in which A is known This method of data ana 
vsis seemed to offer the advantages of simplicity and exactness 
within the limit of the assumption of constant ( in that no 
terms of a series were neglected However, the time-distance 
data must be converted to velocity-distance data which involves 
a hidden pitfall. The velocity obtained by dividing a given dis- 
tance interval by the corresponding time interval is not the ve 
locity at the center of the interval because velocity is not linear 
with distance 
The point within the interval at which the average velocity in 
the interval equals the local velocity must, therefore, be found 
) 


We integrate Eq 2) to obtain 
V = Vie~KCo 3 


where JV, is the initial velocity \ second integration over the 


length and duration of the interval, Y and 7, gives 


T = (1/KCpV,) (e®©2* — 1) = (1/KCp) [(1/V,;) — (1/1 
} 


With this, we can form the average velocity, \, 7, and set it equal 


0.532 


0.528 + 


0.524 - . 


0.520 - — 


n O516 t . 


0.512 t . . 
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Fic 2. One set of data showing the variation of V and 
In [1 + (A/V,)] with x 


to the local velocity at the point »Y where u is to be determined 
For constant deceleration, » would be equal to 1/2 
V = X/T = KCpX/[(1/Vs) — (1/Vi)] = Vie~SOo"X 
Solving for 7, we obtain 

In} [Vi/ Vy — 1] /ln(Vi/ V7) } 


In(Vi/V; 


which has been plotted in Fig. 1. 

To compute a drag coefficient from range data, we would then 
proceed as follows. Choose a set of intervals of equal length, for 
which » is nearly constant, and plot V as a function of x, assuming 
n equals 1/2. (The effect of this assumption for equal intervals 
is primarily to translate the data along x Now, reading the 
velocity ratios from this curve, obtain the exact values of » and 
make the plot of /n V against x to obtain the straight line whose 
slope equals —ACp. An example is given in Fig. 2. Least 
squares fitting of the data to the straight line may be easily ac 
complished when desired. Experience has shown that this 
method offers considerable advantage over those used earlier, par- 
ticularly when the velocity variation in the range is large 

A final point concerns computative accuracy Because /n V 
changes very slowly with V when V is large, we express V in terms 


of a reference velocity and an increment 
V= V+ A= [1 + (4/V,)] V; 
and note that the slope of the line /m [1 + (A/V,)] versus x is the 


same as that of the line /” V versus x 
ln V = In [1 + (A/V, + /n V, 
d(ln V)/dx = d}ln 1+ (A/V, i dx 


Since A/V, is small compared to 1, the logarithm of 1 + (A/V, 
is sensitive to changes in A, and this fact is of some value in further 
improving the ease of calculation 

It has been called to our attention that a somewhat similar drag 
reduction procedure was applied in Germany around 1911 to the 
study of the drag of airships in coasting flight Time rather than 
distance was selected for the independent variable, and, by plot 
ting 1/V as a function of time, a straight line was obtained. It 
can be readily demonstrated that such a presentation will give a 
straight line if Cp is constant. In reference 3, velocity and time 
were the measured quantities, so that the conversion of time-dis- 
tance data and the calculation of the factor m were not required. 
Such a conversion could be made along the lines described herein 
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if the velocity-time coordinates were for any reason preterred 





velocity-distance coordinates 
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On the Buckling Problem of Shallow Spherical; — 
Shells Under External Pressure* 
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— suBJECT of the buckling of shallow spherical shells su it 
jected to uniform normal pressure was recently considered i repo! 
several papers. Two different interpretations of the apparer 
scattering in test data for the critical load g, were given. In pp 
paper by Klein,' it was suggested to attribute this scatter t prox 
initial irregularities of the specimen. On the other hand, it was} aecu 
attempted by Reiss? to interpret the scatter solely on the basis - 
sudden changes in the shape of the normal deflection of the shel lowit 
More precisely, it was observed that there is an oscillating b ‘onl 
havior of the dimensionless buckling load PAw) = de/Ves wher } 
Q. is the buckling pressure of the complete spherical shell fla 
given by the linear theory and yp? is defined by Eq. (2), with shar . 
peaks at certain transition values, u,2 (k = 2, 3, , to be de} the 
fined below. The latter interpretation is supported by a recet defle 
theoretical investigation by Reiss, Greenberg, and Keller r 1] 
purpose of this note is to present some results that differ from} “cla: 
those obtained by previous approximate solutions and that te II 


to support the remarks in reference | 





The differential equations for finite deflections of shallow spher 


ical shells of uniform thickness h, if restricted to axisymmetri If 2 
modes of deformation, can be written in the following dime tion 
sionless form®* § | defle 
LF + w’°G 2x + yFG/x = 0 
{ F 
ae a ae 
LG weF — yF?/2x = 0 | | taine 
T 
where LY = d[x~'d(xV)/dx]/dx, and x, F, and G/x are the di 


mensionless polar angle, angular deflection of the middle surface 


and meridional stress resultant, respectively The geometrica 





parameter uw? whose value is zero for a flat plate, and the load 





° t1o 
parameter y, proportional to the external pressure g (g > 0 in 
ward), are defined in terms of the radius of the middle surface 
° . ‘ tel 
the opening angle 2a, and the elastic constants E and vy as follows 
Kcet 
= [12(1 — v?) ]?/? atatg/4Eht, wv? = [12(1 v?) |? aa?/h (2 ess 
\t the edge, the shell is assumed to be either simply supporte — 
Case S) or clamped (Case C), which can be expressed in terms “ 
Li¢ 


of Fand Gas follows: 


F’(1) + vF(1) = G’'(1) — vG(1) = 0 Case S 





In both cases the conditions of symmetry require 


* The present note is a report on work done under the sponsorship of the | ert 
Office of Naval Research, carried out at the Computation Center of the seem 
Massachusetts Institute of Technology, and supported, in part, by th 
International Business Machines Corporation it 

The writer wishes to express his appreciation to Professor Eric Reissnet valu 


for his helpful advice and guidance during the course of the work 
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FiO) = G(O) = 0 j 
\pproximate solutions of this nonlinear boundary value prob 


were obtained in Case C by Kaplan and Fung* and by Archer 
rhe 


obtained by 


the basis of perturbation methods “exact”’ solutions 


this note have been assuming power 
series expansions for F and G, F = 3,” F,.x", G = 2,°G,x". By 
to suitably 
get 


for a wide range of the parameters y and 


reported 1 


ipplying an iterative technique chosen initial ap 


proximations, it was possible to solutions of a prescribed 


curacy (0.01 per cent 
uw. «Th 


lowing discussion is restricted to the determination of the dimen- 


e results are given in more detail elsewhere,’ and the fol 


sionless buckling load p, as a function of yu? 


the shell deforms similarly to 


For small w?—i.e., 0 < pw? < uw)? 

flat plate with a monotone increasing load deflection curve 
= g( I IW” = normal displacement of the middle surface at 
le apex x = 0 If w? is in the interval w;? < uw? < we, the load 
deflection curve g = g(W) possesses a relative maximum at HW = 
Wand a relative minimum at HW’ = W,, > W \ccording to the 
classical’ buckling criterion, the critical load is given by g. = 
H More generally, we define the buckling load g. as the 
irgest load for which the boundary value problem constituted 


yy Eqs. (1 3), and (4) for fixed w? has more than one solution, 


If w? is in the interval u < pw? < pw”, the shape of the load deflec- 


tion curve g = g( IW) is more complicated. Moreover, the normal 
deflection w(x) changes its ‘‘mode,”’ that is, the maximum of 
does not occur at the apex, as is the cas¢ for pe ~~ 2 
For a limited range of w?—1i.e., for p* 20 in Case § and for 
ue I6 in Case C load deflection curves g = g(IV) were ob 
ined also in the post-buckling region —i.« for deflections 
J I”. and even for loads For uw? = 20 (Case S) and 
u? = 25 (Case C), the results obtained for ulow making a 
irly accurate estimate of g- by extrapolation. However, to de 
rm for larger values of uw? —e.g., w? = 49, Case C— calcu 


lations must be carried out for values of y of the order of magni 


tude 1000, which in the present method of solution would require 


retention of a number of terms in the power series that would 


xceed the capacity of the computer (Whirlwind-1). Dimension 


ess buckling loads p. = QO. = yc/u' obtained for uw? < 20 
Case S) and for uw? < 25 (Case C) are plotted against yp? in Fig. 1 
rhe broken lines indicate the largest values p = g/Q., / p., for 


which calculations have been successful by the present method 


for 25 as lower 


r wu 2), therefore be considered 
Note that 


nin reference 3 are considerably below these bounds 


and which must 


bounds for the actual values p some of the values 


\lthough further computations would be necessary to extend 


the rves p.(u2) bevond pw? = 25, the results show that the os 
latory behavior of p-(u2) with maxima near wp? = ux? (k 2 1) 1s 


as pronounced as is assumed in reference 2. It 





not 
seems somewhat fortuitous that the peaking behavior due to the 


experimental scattering shows the same trend as the theoretical 


given in reference 3, which in turn appear to be, at least 


Values / 
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in part, a result of the failure of the computing process to con 
verge to an expected solution 
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Effect of Dissolution on Oleopneumatic 
Shock-Strut Performance* 


Charles W. Bert 
Principal Mechanical Engineer 


Battelle Memorial Institute, Columbus, Ohio 


ober 14, 1957 


Oc 
ict 


uircraft-landing-gear strut, damping is 


I AN OLEOPNEUMATIC 

provided by forcing oil through an 
a volume of a1 
load-stroke¢ 


theory and 


irifice, and 


rhere | 


springing 1s 


achieved by compressing is been good 


igreement between pneumatic performances pre 


dicted by simple air-compression ictual performance 


under rapid loading, such as experienced in landings and drop 
tests at high sinking speeds.' However, there has been a dis 
slowly 


crepancy between theory and experiment for static and 


is found in gentle landings and 


ipplied loading, such taxiing.’ 


\n analvsis which accounts for this discrepancy is presented 
her¢ 

\ useful relation between the instantaneous absolute pressure 
p and the corresponding density p of a constant mass of air in 
evlinder is the familiar equation 

p/p’ constant l 

where ? p rly tropic-process exponent unit for static load 
ing. Walls! found that the dynamic value of 7 iried from 1.01 
to ] 10 

It is well known that in a hydraulic oil system under pressure, 
i certain portion of any air present in the system becomes dis 


solved in the oil after a sufficient length of time. This phenomenon 
may be visualized physically by thinking of some of the air 
molecules being interspersed between the oil molecules. Be 


cause of the presence of air and oil in physical contact with each 
other in a shock strut, it is reasonable to expect dissolution to take 


pl ice there also 
The 


volume V 


dissolved air volume Vy is a function of pressure Pp, oil 


, temperature, and time of exposure of the two fluids to 
Campbell 


value of Vg, measured at 


each other presented the following expression for the 
saturation standard temperature and 


pressure, 
V KV. 2 

where A a constant depending upon temperature 

ill of the 


In using Eq. (2), it is assumed that the oil contains 


Fairchild Aircraft Division 


Md 


rhis problem was suggested by J. D. Stites 


Fairchild Engine and Airplane Corp Hagerstown 
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Strut Position X, Measured From Fully Extended Position, in 


Fic. 1. Static load-stroke curve for a typical shock absorber. 


air which can be dissolved at any particular pressure. Thus, the 
time lag between air compression and dissolution, or between air 
expansion and effervescence of the dissolved air, is assumed to be 
negligible. This is probably not valid for rapid loading, as shown 
by data obtained by Walls on air compression in which the total 
loading time was less than 0.25 sec. However, there are indica 
tions that dissolution may be taking place in some gentle drop 
tests,? and data to be presented later verify that this effect takes 
place in static loading 

The total mass m of air, including both gaseous and dissolved 


phases, in the cylinder is given by 
m= p Vv. + p V = px V. + AK V.p 3) 


where subscript s denotes quantities at standard conditions 
When the cylinder is under a pressure p, the density p of the air 


in the gaseous phase is 
p=i(m pKVip) A/V 1 


where | = air-space volume corresponding to pressure p 
\fter substituting the right-hand side of Eq. (3) for m in Eq. 
#) and applying Eq. (1) with » = 1 to the air in the gaseous 
phase at pressures p and p,, the following result is obtained: 
pV = pslV. — KV.(p — ps)] 5) 
The pressure p;, and air-space volume V,, at the fully extended 
strut position may be substituted into Eq. (5) to vield 


Ve = Vee (pre/Ps) + KVolpre — Ps) 6) 

Substitution of V. from Eq. (6) into Eq. (5) gives 
Pp = [(Vre + KVops)/(V + K Vops)] pre 7 
Now, 41 = piston area, and x = any strut position, measured 
from fully extended position. Substitution of (V;, —- Ax) for V 


and strut loads ? and P,;, for pressures p and p;, in Eq. (7) results 


In 
P = [(Vie + KVips)/(Vie — Ax + KVip.)| Pye 8) 


\lthough Campbell gave 0.0073 psi”! for A for saturation of 
air in oil, he gave no information on the time required for satu 
ration or the effect of temperature. Therefore, the static value of 
K was determined during compression loading of a strut in a 
universal testing machine. Design data included: A = 10.31 
sq.in., Vo = 167.3 cu.in., V;- = 121.5 cu.in., P;. = 2,940 Ibs 
Test results, plotted as the reciprocal of the load versus the 
stroke, are given by the encircled points in Fig. 1. It is to be 
noted that they fall very close to a straight line which differs 
significantly from simple air-compression theory. The straight 
line fitted through the data points was used in conjunction with 
the design parameters given above to determine the effective 
value of A in this test, 0.0038 psi 


1 


Although considerable dis- 
solution took place, the amount was less than Campbell's satu 
ration value. This difference may have been due to a difference in 
the temperature of the tests 


REFERENCES 


Walls, James H., Investigation of the Air-Compression Proce During 
Drop Tests of an Oleo-Pneumattc Landing Gear, NACA TN 2477, 1951 


AL SCIENCES FEBRUARY, 1958 


2 German Landing Gear Design and Testing, Aircraft Eng Vol. 2 
134-136, 150, May, 1948 

3 Campbell, James E., Investigation of the Fundamental Charact ti 
High Performance Hydraulic Systems, AF TR 5997, Air Materiel Commar 


Wright-Patterson AFB, June, 1950 


Analytical Results of Certain Nonlinear 
Flutter Problems 


S. F. Shen* and C. C. Hsu** 

Department of Aeronautical Engineering, 
University of Maryland, College Park, Md 

October 9, 1957 


F_NHE RECENT WORKS of Woolston et al.':* reported on ¢ 
r behaviors of wing-control surface flutter systems with sever 
kinds of structural nonlinearity as solved by an analogue cor 
puter. A number of interesting phenomena were noted such 
the occurences of “limited amplitude flutter,’’ ‘violent utter 
and “‘divergent flutter’? depending upon the system and t 
implitude of the motions. It is felt that apparently if the no 
linear flutter problem can be handled analytically, even in 
approximate manner, much light will be shed on the roles of t 
various types of nonlinearity in producing the exhibited nor 
linear effects 

We have applied the well-known theory of Kryloff and Bog 
liuboff* in nonlinear mechanics to some of the flutter problen 
treated in references 1 and 2. Essentially, at flutter, a stead 


oscillation is assumed. The nonlinear element will generat 


higher harmonics in addition to the component of the funda 


mental frequency. By balancing the fundamental harmoni 
only, the problem is reduced to an equivalent linear one wit 
certain elements depending on the amplitudes of motion as paran 
eters. This is nothing but the simplified version of the ‘‘princ 
ple of harmonic linearization” of Kryloff and Bogoliuboff 

By such procedure, all the structural nonlinearities (excey 
hysteresis) of references 1 and 2 are equivalent to linear stif 
nesses which, however, vary with the amplitudes of motior 
The flutter speed and frequency, therefore, can be directly deter 
mined on the basis of the usual linear analysis. The resul 
naturally depends upon the amplitudes of motion. In the hys 
teresis case, the equivalent linear stiffness is complex—.e., a cet 
tain amount of structural damping occurs. Again, both the re 
and imaginary parts vary with the amplitude 
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Fic. 1. Bending-torsion flutter with free play and preload | 
the torsional stiffness. 








FIG 


el Commap 


ir 


ted on t 
vith sever 
logue cor 
ed such 
it Nutter 
n and t 
if the nor 
‘ven in 

oles of th 


bited nor 


and Bog 
problen 
a stead 
generat 
he funda 
harmoni 
one wit 
as paran 
e *‘princi 
ff 
Ss Cexcey 
lear stiff 
f motior 
tly deter 
he resul 
the hys 
e., a Cer 


h the re 


READERS 








37 














1-0 If 1-2 


V/VuIN 


Wing-control surface flutter with free play and pre 
load in the aileron stiffness 


\ few of the cases in references | and 2 are worked out by the 
present method in reference 4. The case of wing bending-torsion 
flutter with free play and preload in the torsional stiffness is 
shown as Fig. 1. The exact value of 6, was not given in refer- 
ences 1 and 2. We 
ér is kept at 0.5° as in the original references 
metric stiffness property with regard to the direction of motion 


have computed for 6, = 0.2° and 0.5°; 


The nonsym 


causes a steady deviation, effectively moving the neutral posi- 








load 


tion during oscillation a small distance toward the ‘“‘soft’”’ side. 
It does not affect the flutter determination 
According to Fig. 1, our analysis indicates that for V/ Viin < 1.0 


thesystem is stable for small amplitudes but may assume very large 
amplitudes if the ratio is increased toward unity. Such behavior 
seems to agree with the description of references 1 and 2. How- 
ever, the flutter boundary of Woolston et al. (Fig. 1 of reference 2) 
as superposed on Fig. | is seen to be rather different in character 
We are quite certain that, as the amplitude increases much 
larger than the value 67, the system certainly should behave again 
as the linear one in the manner predicted by our analysis 

In Fig. 2 we show our results for the case of bending-torsion 
control surface flutter with free-play and preload in the control 
)-. 


stiffness. The ratio of 6,/67 again was not given in reference 2; 


our curve is for 6,/67 = 0.5. The region to the left of the bound- 
ary is stable, to the right, unstable. For 1.2 > V/ Viin>1.0, 


we see that the system is unstable for small amplitudes but will 
stabilize at an amplitude of 8/57 < 2; vet, by further increasing 
Beyond V/Viin = 


This prediction 


the amplitude, divergent motion occurs 
1.2, the system is unstable for all amplitudes 
also agrees well qualitatively with the result of reference 2 (Fig 
6 \gain the nature of the flutter boundary is quite different 
To us, ours seems to be more convincing. In fact, it appears 
questionable whether the scale for 8)/é7 of Fig. 6, reference 2, 
might not be erroneous. For 80/67 of the order of 30, reference 2 
indicated a much stronger nonlinear effect than our prediction 
for 8/67 of the order of 5. To us, such a persistence of the non 
linear effect is rather difficult to accept 

In reference 4, the case (A) of Fig 
While qualitatively in agreement with our analysis, 


the numbers given by reference 2 must not 


5, reference 2, was also 
examined 
in this case at least, 


be accepted without reservation Even for 8)/67 = 0, the value 


of V/Viin was given as approximately 1.75, in contrast to the 
value 1.46 directly obtainable from reference 5 
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Stability Equations for Conical Shells * 


Joseph Kempner 


Professor of Aeronautical Engineering, 
Polytechnic Institute of Brooklyn, Brooklyn, N.Y 


October 9, 1957 


O* Pp. 47 OF REFERENCE | Vlasov presents equations suitable 
for stability analyses of arbitrary thin-walled shells of con 
For conical shells of circular cross section, these 


stant thickness 


equationst become 


DV *x F..,/x) cota Ne 
NO 9[(w,99/X? sin? a) 4 
2.N°.9(w.9/X Sin a + Vu + Ow y X SIN a 0 l 
ViF + (w,,,/x)Eh cot a 0 2 


Here w is taken as the component of the buckling displacement 
normal to the median surface and directed toward the anis of the 
cone whose apex angle is 2a and wall thickness is h; x is a co 
ordinate measured from the apex along a generator of the median 
measured in planes 


surface; and #@ is an angular coordinat 


normal to the axis (x, 3, and the inner normal form a righ 





handed coordinate system). N®,, N®s, and N°,» are t 
brane stress resultants prior to buckling (their sign convention 
corresponds to that of elasticity theory Y and 9, respectively, 
are distributed surface loads fixed in the axial and circumferen 
tial directions (together with applied radial and edge loads which 
do not appear explicitly in Eqs. (1) and (2), they are considered 


in determining solutions of the membrane equations mentioned 


previously D = (1/12)Eh?/(1 v?), where E is Young’s 
modulus and vy Poisson’s ratio. The operator V‘ is defined as 
vi = V2V2, where V*( = (1/x)|x( + (1/x? sin? a ry) 


and the stress function F is related in the following manner to the 


increments in the stress resultants arising from the buckling 
displacements: 
NV, = (F.95/x? sin? a) + (F,,/n Vy = I 
V9 = f sin a 
It can be readily shown that 
1/x) (xtViF = H(F 3 
in which //7 is an operator defined by 
H() = (1/x)V¥* + (8/x? + (24/2 
SX + S A s1 t 

Eqs. (2) and (3) yield 

H(F..:) = 1 /x \ Eh cot a j 


ind, therefore, from Eqs. (1) and 


* This work was performed under a consulting contrac 


Electric Company, Missile and Ordnance Systems Department, Philade 


phia, Pa 

+ After typographical errors on p. 47 of reference 1 are removed, the coeffi 
cients of the first quadratic form are taken as A land B sin a, and 
it is assumed that the stress system prior to buckling is governed by the well 
known membrane equations for conical shell 
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H(xV tw — (x/D) | Nw + N%y[(w,99/x? sin? a) + (w,,/x)] + 2N%,9(w,59/x sin a), — Xw,, — O(w.9/x sin a \) + 
Eh D) cot? a hy ye 0 


This eighth-order equation in w can be considered to replace Eq. (1) which is coupled in w and F. An equation similar to Eq 
was derived by Seide in reference 2, but his equation differs from the present one in several terms in 7. It is believed that this is simy 
due to an error made in the determination of a relation in reference 2 similar to Eq. (8 

At times, it is convenient to have the stability equations in terms of the buckling displacements. From pp. 32 and 47 of referenc 
it can be shown that the following two equations represent conditions for equilibrium of the membrane stress resultants arising duri; 


buckling 





9 3-1 l l y /X — 12/(1 + 1 ) [2 t/Xx vj; sina — 
l— vp 1+» U 39/X sin a 7 2 ; 8 vu, = W/X COS @ 
{ ] >— 7 | u Hiy/ Xk — 2 u aa /X sin a — 
l— 4 1+ v)] } [x3(o/x vi sina + [2/(1 i »/x) cota (7 
in which uw and 7, respectively, are median-surface buckling dis 
placements in the x and J directions. Eq. (1) represents a third 
tion in the buckling displacements sin Taide, Pad, A Donntlype 7% apn # 
‘ 0 > uc . < spiace “: 2» > ce 
ce : iinet. : Buckling of Thin Conical Shells, Guided Missile Research Division, The Ran 
F.. = Ny = [Eh/(1 y2)1[(1/x sin a)(vys 4 Wooldridge Corp., Report No. GM-TM-103, July 13, 1956 Als 
: : ' Appl. Mech. Paper No. 57-APM-42) 
Wu sin a@ w cos a) + vu : 
Hoff, N. J., Thin Circular ¢ tcal She Und irbitra Loa A 
in accordance with reference 1 Mech., Vol. 22, No. 4, pp. 557-562, December, 1955 | 
° | 
If it is assumed that 
| 
u,z>>(u/x) = (tan a/r)u, v,2>>(v/x) = (tan a/r)v (8 + 
where ¢ is that principal radius of curvature of the median surface 
which lies in a plane normal to a generator, Eqs. (6) and (7 . 
a ' lads is bei Some Remarks on the Foundations of 
can be readily transformed into more convenient expressions . . . 
a 5 sige nae I oi Piston Theory and Supersonic-Hypersonic 
le assumptions made in Eqs. (8) are certainly reasonable for H H H * 
I Similarity 
thin-walled, truncated cones of small apex angle (see reference 3 
furthermore, as @ — (7/2) the cone approaches a flat circular is 
| only Eq. (1) (which 1 cag et | F and ish 
Nate and only u¢ (which becomes independent of / and, r 
. : / ot Aerodynamics Division, Weapons Research Establishment 
hence, of « and 7) is required Salisbury, South Australia 
If Eqs. (8) are applied, Eqs. (6) and (7) reduce to October 14, 1957 
UV. = 2/(1 + v)|(x sina u , (1 v)/(1+pv)] X 
(u,99/X sina) + [2/(1 + v vw, (w/x)| cosa (9 I {** s!' SHOWED THAT when the velocity changes in the fre 
me : stream direction are small in comparison with those in t 
Urs = 2/(1 + v)| (v7. 99/X sin a é : a : : 
: . ‘ transverse directions the equations of motion of steady compres 
(1 — v)/(1 + v)| (x sin at + , . : : 
sible flow approximate those of unsteady flow having one les 
2/(1 + vp (ww 9/X) cot a (10 5 “ . ‘ : : 
dimension. Considering the hypersonic flow past a two-dime: 
With a procedure analogous to that given in reference 5, these sional airfoil, the above condition is nearly satisfied and the fl 
equations transform into pattern becomes analogous with that produced by a piston in tl 
a Co me ae —_— 4 distance-time plane. The nondimensional piston velocity cor 
Hy: cos a }(2 v) X (XW,79),, Sin? a@ Ww 9o9 list t l rl l l t locit 
be eh htt v)|xw.r9 sin? a} (11 responding to a given flow deflection angle (6) and free-streai 
Hu) = cosa sin a {oxlx(xw . sin? a cai Mach Number (.1/) is 1/7 tan 6 (or \/J7) which is the similarit 
; , parameter proposed by Hayes and, prior to him, by Tsien.’ 
Vw sin? 2/(1 y) |W ve (12 : . - 44 
, . : wey However, the approximation made by Hayes can be avoided in | 
in which //; is an operator defined by some cases since the equations of motion of two-dimension 
Hu ox { x(x ) ).. sin! steady flow reduce to those of one-dimensional unsteady flow i 
a v) ‘ > aT ° . P . . 
the velocity component in amy direction remains constat 
7 a2 n . " 
aN IN roe SI” a III throughout the flow field When such a direction exists, t 
(see reference 3). Except for the existence of the last term equivalent piston velocity is in a direction normal to it and 
within the braces in Eq. (11) and the last two terms within the magnitude such that the body boundary condition is satisfi 
braces of Eq. (12), these equations are the same as the corre- The nondimensional piston velocity which is the ratio of the pis 





ton velocity to the speed of sound in the undisturbed gas the 


sponding equations of reference 3 


° A - ° oe ‘rve j ; -} te 7 } c , 
Expressions relating membrane, bending, twisting, and trans serves as a similarity parameter for the steady flow 


: » case of . ic o-di nsi ows (obli shock 
verse shear stress resultants to the displacements a, 7, and w can In the case of the basic two-dimensional flow blique shioc 


‘ - oe » . oye Ss ) ] ' S S ich ce sis SS i 
readily be obtained from pp. 27, 28, and 32 of reference | These waves and Prandtl-Meyer expansions) which consist essentially 


relations (together with Eqs. (8) if Eqs. (11) and (12) are used two regions of uniform flow separated by a transition region, 


must be considered when boundary conditions are enforced direction in which the velocity component is unchanged does exis 


Eq. (11) which relates 7 and w, and Eq. (12) relating u and w outside the transition region and in the oblique shock case it 1s 


together with Eq. (5), represent equations alternative to those of course, the direction of the shock front It should therefor 


formulated previously in terms of w and F Of course Eqs be possible to find for these flows an exact similarity paraimete! 


‘ me . : ic for lo Ss . s high s rsonic speeds d 
11) and (12) depend upon the suitability of the assumptions which holds for low as well as high supersonic speeds and I 
large as well as small defiection angles 


represented by Eqs. (8), while Eq. (5) remains as originally 
derived. Inasmuch as experience with buckling problems of Considering first oblique shock waves (Fig. 1) with the axes 


circular cylinders indicates that « and v play a minor role com- as shown, the flow is analogous to the flow in the (x, ¢) pla 


pared to that of w, it is believed that little accuracy is sacrificed which would be produced by a piston moving with uniform velo 
if Eqs. (11) and (12) are employed ity along the line 48. The boundary condition is satisfied | 
the ratio of the piston velocity (U’) to the free-stream flow veloc 
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ity is made equal to the ratio of the length of AB to that of 
)4 in Fig. 1—1.e 

U/q = AB/OA = (U/a)/M 
Therefore 
U/a = M sin 6/cos (a — 6 l 


where u is the speed of sound in the undisturbed gas and a is the 
tion in which the flow velocity remains constant 

For Prandtl-Meyer expansions (Fig if the convention is 
| that angles are positive when measured counterclockwise 


he free-stream direction, it is found that the nondimensional 





piston velocity is again given by Eq. 1 
Hence for both cases the exact similarity law becomes 


pb f;{.M sin 6/cos (a — 6 


Unfortunately, a cannot be simply expressed as an explicit func- 
tion of V7 and 6 for either of the flows. However, the angle of 
inclination of an oblique shock wave to the body surface is very 

irl\ equ il to the free-stream Mach angle be so that the re- 
placement of a — 6 by wisa good approximation for oblique shock 
vaves. Fortunately, it is also a reasonable approximation for 
Prandtl-Meyer expansions. The similarity parameter then be- 


mes, 


VW? sin 6/Y AJ? — 1 


py = fo( M sin 6/cos p) = 


rin terms of thickness-chord ratio and Mach Number, 
p/p, = fo( MP?7/V M?-—- 1YV 7? + 1 


It has been shown’ that the replacement of \J by 1/?7/Y 1/7? -— 1 
in the conventional hypersonic parameter .J/ tan 6 increases the 
wecuracy at low supersonic Mach Numbers and the use of sin 6 
instead of tan 6 greatly increases the accuracy at large deflection 
ingles 

Piston theory may also be applied to axisymmetric flows and 
the flow over a cone becomes analogous with the flow produced 
by a uniformly expanding cylindrical piston. However, the 
shock wave produced by a cone does not separate two regions of 

mstant flow velocity and direction, and there is no direction in 
that an exact 


which the velocity component is unchanged, s« 
similarity parameter does not follow. The direction of the sur 
face of the cone is found to be the closest simple approximation 
to such a direction, and, when the expanding cylinder analogy is 


ipplied in this direction, the similarity parameter J sin 6 is ob 


tained 
Van Dyke? has proposed the similarity parameter C, = r*f; 
rV iM? — 1) as a combined supersonic-hypersonic rule He 


ioted that p/py = f VW2r/VY M?* — 1) had also been proposed as 
icombined similarity parameter and states that it is restricted to 
planar systems Further, he remarked that hypersonic and 
ombined similarity rules apply to rather thick bodies if 7 is 
taken to be the sine of the maximum angle of inclination of the 
wdy surface. This note gives a simple explanation in terms of 


piston theory of both these observations 
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SUMMARY 


he variational principle of strain energy used in order to determine 





when the turbine is at rest, the direction of the eigenvector in which a group 
of blades with constant cross section and twisted along their length vibrate 
in the fundamental bending mode when they are tied together with one 
lashing wire Investigation shows that in the case of one connection the 
eigenvector direction almost coincide with the direction of the ax of 


maximum moment of inertia 
DISCUSSION 


N IMPORTANT PART concerning the blading problems in 
i \ turbodesign belongs to the determination of the angle in 
which vibration of a group of blades occurs when the turbine 


is either at rest or in motion. Such an angle is called the vibra 


tion angle, 8, and is measured with respect to the major axis of 
inertia. In the literature, the investigation of this problem is 
Mather 


interpretation of this problem implies two separate « 





due to the work of Geiger' and Stanisi 








systems, namely, Ogx and Okx, with O at the intersection of 


principal axis. From this, ind k&-k are axes of maximum an 


minimum moment of inertia, respectively; y-y, direction of the 


blade perimeter; direction of the turbine axis, and x-x, 
direction of the blade axis; a, angle of twist; and #, angle 
between plane of load and g-g axis. For all other symbols see 


Figs. 1, 2, and 3 
Assume that the actual angle of vibration of blade 
coincides with that generated when the blade carries 


} ] 
ya 


load P in addition to a single load Q acting at the | 


uniform 


the lashing wire is hard-soldered to the blade hen, the strain 


energy of the system under consideration is 









TOURNAL OF THE 

















where 


A, = (EI,/2 


aL 
ky. rr(x, 3)dx 4 } ko 2x°A(Xx, 0 is| (3 


are the strain energies of the blade corresponding to the bending 
vd), k(x, 3) 
to each axis, index 1, 


domains 0 < x < hh, 1; << « < Las shown in Fig. 1 
Ag = Mek, 


ii 


are the corresponding 
2 refer to the 


And 


about g-g and k-k axes, g(x, 


deflections with respect 


r(1h, dv) (4 


is strain energy introduced by the bending moment transferred 


to each blade by the lashing wires on both sides. Note that 


QO/P 


AERONAUTICAL 





u = 
Then, using ordinary flexural theory one obtains for 0 < x </ 
(x, 0) = (x?/2EI,)P\(L/12) [6 1(x/L) + (x/L)?] 
ull: — (x/3)]} sin d (5 
k(x, 3) = (x?/2EI,)(P}(L/12)[(6 — 4(x/L) + (x/L)?] 4 
ull, — (x/3)]{ cos 3 2M) (6 
and for]; << x < L: 
g(x, 0) = (x?/2EI,)P}(L/12)(6 4(x/L) + (x/L)? 
u(l;/x)?[x — (41/3)]f sind (7 
ko(x, 9) = (x?/2EI,)(P}(L/12) (6 4(x/L) + (x/L)?] 4 
u(l,/x)2[x — (1,/3)]} cosd — 2Mpz) (8 
TABLE 


Geometrical Characteristics of Turboblades for Driving the 


i l t hag F Fz 

in in.) (in.) (in.) (in.? in.” 
J-33 $00 2.66 Li 0.80 0.279 0.020 10 
J-35 3.50 2.33 0.58 0.29 0.197 0.020 10 
J-40 3 §2 2 Ge 0.62 0.22 0.138 0.020 10 


SCIENCES FEBRUARY, 1958 


From Figs. 2 and 3 it follows that 





ki, 


3) = 2 COS a 
Evidently 
= jf,/* LDLtEI 
fo = (1/2)t sin ag L,, J 
From Eqs. (6), (9), (10), and (11) it follows that 
We = Al([p/I;,) cos 3 cos a (1p/I,) sin 3 sina 9 


where 
3LP{1L + (1/L)(u 1) + (1/3)(),/L)? 
12(/ 2/1 


} ] #2 > 
(i yt COS a 





Hence, Eq. (1) can now be evaluated. It follows that 


A(J)) = R, sin? 3d + Rs cos? 3 R; sin 26 
where 
R, = (L*P?/8EI,), (1/5) + (4/L)?[2 $/3)(1,/L 
1/3)(1/L)? | (4/3)(1/L 33} 
R. = (L'P?/8EI,,.)) (1/5) + (1,/L)?[2 (4/3 LE) 
(1/3)(1,/L)?]u + (4/3)(/L 374 (PAL/2EI,) X 
(p/h il + h/ Lu l 1/3)(1,/L)? 
(2/3)L{1 + 8u(h/L)2|} cosa (I 
R; = (PAL/2EI,)([e/1,),{(1/3)L[1 3u(1,/L)? 
(1/2), [1 + (11/L)(u l (1/3)(1,/L)?]} sin a 7 
Note re. resulting de flection of the blade, then 
f(L, 3) = gL, 3) + RAL, V 18 
or f-(L, 3) = U, sin? d + Ls cos? i + L3 sin 2d 19 
where 
U, = (L3Pe¢/2EI,)* + (L7A/EI;,)*(1p/1,)? sin? a 2 
Us. = (L*Pe¢/2EI,)? + (L?A/EI,)*( 13/1)? cos? a 
LiP Ae (E] 2 (Ip I COS a PA | 
U; = (1/2), [(L>PAe/(EI,,)?](1e/I,) sin a 
(122A /EI;-)?(1g?/Ij1y) sin 2a} 22 
e = (1/4) + w(l,/L)? [1 (1/3)(4/L | 


Now the strain energy per unit deflection can be obtained 
The necessary condition from which the deflection angle 0 can 
be obtained follows from the variational principle of the strain 


energy per unit deflection in the system—i.e., 

A(d)/P(L, 3) = Minimum 23 | 
or A gd) PCL, 3) — fe(L, IACI) = O 24 
Eqs. (14), (19), and (24) lead to 








LASHING WIRE 








Curved part of lashing wire corresponding to blade 


Fic. 3 


l 
Compressors J-33, J-35, and J-40 


i. I, Ip a v p 
(in. 4) in.* in.! (rad.) (rad 
6 X 87210 10-§ K 482 1078 & 30 95°20) 0.00198 0.00 
6 X 28230 10-8 & 369 10-6 X 30 x) 0.00442 0.00 
6x 1f 10-8 & 162 10-* X& 30 $2°35’ 0.00825 0.00 
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’r, Sin 20 + (72 + 7 








= U.R U| Re, re = 2UR U:R,), 7; = 20U2R 
R From Eq. (25) it follows that 
} $2\(r2 + ors) VW ri? + orets + ri(t r;)| A sa 
J sin! * (26 
' 47:7 + (72 + 13)? 
But at the end of the blade, the following relation exists: 
f(L, 3) = gAL, 3) sin B + R(L, 3) cos B (27 
Then, the vibration angle is given by 
B = sin v(L, 3)/f(L, 3 (28 
Table 1 shows the results obtained for turboblades for driving 
e compressors J-33, J-35, and J-40, at the laboratory of the 
School of Aeronautical Engineering at Purdue University The 
results show that in the case of one connection the eigenvector 
direction almost coincides with the direction of the axis of 
maximum moment of inertia 
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Concerning the Rigidity of Elastic Structures 
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I IS THE PURPOSE of this note to state and examine some limi 
tations of a theorem often used intuitively by the stress and 


vibration analyst. This theorem can be expressed as follows 


\dding material to any traction-free subsurface of a complex 


elastic structure which is acted upon by a system of external 


forces results in a decrease in the deflection under each load, and 


vice versa. This theorem is often used intuitively when the 


design engineer ‘beefs up’’ a structure to make it more rigid 

Biezeno and Grammel! have introduced the above as an 
iuxiliary principle necessary for the deduction of Saint-Venant’s 
However, instead of the actual displacement, they 


The 


theoretical development leading to inequality (5) of this note is 


Principle 
use an “average’’ displacement which is not clearly defined 


different from their physical reasoning and in addition deals ex 
plicitly with the actual deformation and stress fields 

as shown in 
The 


* and e€ 


Let us consider a general structure (structure 1 


Fig. 1, excluding region A having external tractions 7 


associated internal stresses and strains are given by o 


Further let us consider a second structure (structure 2) which is 


identical to that indicated above but in addition has material 


added to some traction-free surface of structure 1 (say region A 
rhe tractions, stresses, and strains in this second structure are 


* 


designated by 7;*, oi;, and €;, respectively 


Let us consider fictitious stress and strain fields o;;** and e€;;* 
in structure 2, made up as follows: 
;;** = oj;* fin the region circumscribed by structure 1} ' 
aa ( 
0 = 0 fin region A{ ; 


¢;;"" = O ) { 


** 


It can be easily shown that o is in equilibrium throughout 


the body, even across the traction-free subsurface S. «;;** 
satisfies Hooke’s law everywhere but will in general violate the 


conditions of compatibility 


RS’ 
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Subject to the above conditions there results,” 


. 

| Gi;" €ij" av — / a, ,€,)di > i ) 
> ) 

J - J 2 


where the integration is over structure 2 


This is valid even for the discontinuous stress and strain dis 


tribution @ and ¢,;**, the only requirement on continuity 


being that equilibrium is satisfied over the entire volume 
Introducing the stress and strain distributions--Eqs. (1) and 


2)—inequality (3) becomes, 


This inequality furnishes a comparison of the total internal ener 
gies of structures 1 and 2 
the actual 


thes 


Since the stresses and strains in inequality (4) are 
values associated with the corresponding structures, satisfy 


equilibrium and compatibility Hence, applying the theorem of 


virtual work there results, 


f Ti *(uz* — udS > 0 5 


where u;* and wu; are the associated generalized displacements 
in the direction of the external tractions 7 

The following general results can be deduced directly from in 
without further calculation 


the 


equality (5 
1. For 


load, adding material decreases the associated deformation under 


system with a single concentrated generalized 


that load, and vice versa. Stated differently, the system be 
comes stiffer when material is added 

2. For the system with a distributed external traction that is 
either zero or constant everywhere else, the average deformation 


over the loaded area in the direction of the traction Ss uid S/S 
decreases with an increase of material 

The results given in 1. above cannot be generalized to hold 
for an arbitrary variable traction or combination of concentrated 
generali 


for 


As a result of this conclusion, each complex 


tractions. To illustrate the impossibility of such a 


zation, examples and counterexamples will be shown some 
simple structures 
loading or class of loading will have to be investigated as a special 


case for the particular class of structures under consideration 


Fig. 2, illustrates a simple support beam with all external 
loads in the same direction. It is apparent that increasing ma 
terial on a traction-free surface will decrease the deformation 


under each load 


Fig. 3 illustrates a clamped beam with a slot parallel to the 


neutral axis. Increasing the slot length L is equivalent to re 


moving material, and, decreasing L is equivalent to adding ma 


terial. For the case of P; > P», there exists a slot length 1 such 
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FIG. 4b 


that increasing L will decrease 62, and conversely, decreasing L 
will increase 62. Thus even though removal of material will in 
crease the total work it will not increase all deflections in this 
case 

Fig. 4a illustrates a structure that furnishes both an example 
and counterexample. For the proper ratio of P;/P», and rigidi- 
ties on each side of the right hand support, the deflections can 
both be made positive as shown in Fig. 4b. Increasing the 
rigidity in region B has the effect of decreasing 6. without affect- 
ing 6;. Increasing the rigidity in region A has the effect of in 


creasing 62 while decreasing 6, 
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The Effect of Internal Pressure on the Buckling 
Stress of Thin-Walled Circular Cylinders 
Under Combined Axial Compression and 
Torsion? 
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October 17, 1957 


| gress AIRCRAFT and missiles often incorporate pressurized 
circular, cylindrical, monocoque fuel or instrument com- 
partments in their fuselage structures. Even though these com 
partments may not be pressurized for structural purposes, svs 
tems requirements often dictate the use of internal pressure. If 
internal pressure is present, it behooves the structures engineer 
to account for the increased stability of the shell (and the sub 
sequent increase in the strength-weight ratio) caused by pres- 
surization. Although previous investigations have been made of 
pressurized cylinders loaded in axial compression alone!:?* or in 
torsion alone,*® no information has been published on the sta- 
bility of pressurized cylinders under combined axial compression 
and torsion 

For circular cylinders which are not internally pressurized, an 
empirical interaction curve for combined compression and tor 
sion has been suggested® in the form 


R. + R? = 1 (1) 

where 
R, = critical compressive stress under combined load divided 
by critical compressive stress under compression alone 
R, = critical shear stress under combined load divided by 


t 
critical shear stress under torsion alone 


t The authors gratefully acknowledge the assistance of the personnel of 
the Structural Test Laboratory of the Missile Development Division who 
conducted the test program 

* Senior Structures Engineer 

** Group Leader, Structures Data and Methods, 
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The theoretically derived interaction curve obtained by so} 


ing Donnell’s equation is almost a straight line’ and has bee 


found to be too conservative when compared with test result 
As a consequence, the empirical curve is generally accepted pr 
vided the values used for pure compression and pure torsion j 
the stress ratios are determined from empirically defined desig 
curves or from tests 

For pressurized cylinders under combined compression ar 
torsion, a theoretical interaction curve cannot be derived fro; 
Donnell’s equation because small deflection theory does not jy 


dicate the stabilizing effect of internal pressure on the compres 


sive buckling stress. To further complicate the problem, no e, 
perimental data for pressurized cylinders under combined load 


have been available in the literature to date. It is the purpos 


of the present note to show that the shape of the empirically d 
fined interaction design curve given by Eq. (1) is valid for pre 
surized circular cylinders under combined compression and tor 
sion 

To obtain test data on the buckling of pressurized cylinder 
under combined loads, a series of tests of circular cylinders wa 
performed in the Structural Test Laboratory of the Missile Dy 


velopment Division \ total of 28 pressurized cylinders wer 


tested, 10 in axial compression, 6 in torsion, and 12 in combine 
axial compression and torsion. The test cylinders were fabri 


cated from 0.0032-, and 0.0087-in.-thick, half-hard, 18-8 stainless 
steel. The cylinders were prepared with a length of 21.5 in., a 
radius of 8.75 in., and had r/t ratios of approximately 1,000 and 


2,730. The cylinders were internally pressurized to 8.4 psi 


Most of the cylinders were fabricated by wrapping the sheet and 


adhesive bonding a single 3/8 in. wide longitudinal seam wit! 


Epon VI adhesive. Some of the later cylinders were fabricated 


by seam welding which proved to be a faster and simpler tech 
nique. 

The test cylinders were clamped by straps to the 3-in.-thick 
circular aluminum head plates of the test jig (see Fig. 1 of refer 
ence 3). The straps prevented slippage and leakage between the 


heads and the cylinder. Internal pressure was pneumatically 





Fic. 1. Typical buckle pattern for pressurized cylinders under 
combined loads 
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COMPRESSION STRESS RATIO, Re 


Interaction curve for pressurized cylinders under com 


Fic. 2 
bined compression and torsion 


supplied, but, for reasons of safety the cylinders were partially 
filled with water. A mercury manometer was used to measure 
the internal pressure. The axial compressive load was applied 
by a hydraulic strut and measured with an SR-4 load cell. The 
torsion load was applied by a hydraulic strut which was attached 
to the lower head by two cables. The torsion load was also meas- 
ured by an SR-4 load cell 

The buckling loads of the test specimens were determined 
visually, which was a relatively simple task because the buckles 
Occasionally, buckles appeared to 


ften snapped into position. 
The order of 


grow out of relatively large initial irregularities. 
ipplication of the compression and torsion loads did not affect 
the magnitude of the critical stresses. A photograph of a typical 
buckle pattern in one of the specimens subjected to a combined 
load is shown in Fig. 1 

The test data are plotted in Fig. 2 in terms of the nondimen- 
sional stress ratios R, and R,. The stress ratios are defined as in 
Eq. 1, but the critical stresses are those obtained with internal 
pressure present. In order to plot the combined loads data, the 
iwerage of the experimental buckling stresses for the cylinders 
loaded in compression alone and those loaded in torsion alone 
have been used as the denominators of the stress ratios. Shown 
in Fig. 2 is the interaction relationship defined by Eq. (1) which 
has been suggested for unpressurized cylinders. It can be ob 
served that this interaction curve is a suitable lower bound to the 
test data. Based on these tests, it is concluded that the same 
shape design curve adequately describes the behavior of both 
pressurized and unpressurized circular cylinders under combined 
For the design of pressurized cylinders, 


should be used with the critical buck 


compression and torsion 
the stress ratios in Eq. (1 
ling stress for pressurized cylinders under compression alone or 
under torsion alone determined from empirical curves 

It should be that the above conclusions are 
based on tests of cylinders at only two values of the pressure 
It is conjectured that the interaction 
long pressurized cylinder; 


remembered 


parameter, (p/E) (1/t)* 


same for any 


curve would be the 
however, additional tests should be performed to verify this 


issumption. 
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SYMBOLS 


step height 


bubble length 


n number of 6.* 

R;.* U s5s* 

Ua free-stream velocity at reattachment 

ls free-stream velocity at sepa: ation 

Xa distance from leading edge to reattachment 
Xs distance from leading edge to separation 

a angle of attack of airfoil 

ba* displacement thickness at reattachment 
6,* = displacement thickness at separation 


kinematic viscosity 
DISCUSSION 


oo THE ADVENT Of airfoil shapes having thin leading edges, 


attention has been focused on the attendant leading-edge 
separation bubbles since they may greatly affect the downstream 
boundary layer as well as the airfoil lift characteristics 

Owen and Klanfer,' based on the results of experimental in 
vestigations, presented a correlation of nondimensional bubble 
lengths with a Reynolds number based on the displacement 
thickness and the free-stream velocity 
They observed, from this correlation, that the bubble lengths 
divided into two general groups. The short “turbulent” type 


500 while the longer, ‘‘laminar”’ 


at the separation point 


bubbles occurred for R5,« 
type occurred for R5,* << 500 

Due to the complex nature of these separation bubbles, only a 
limited number of incomplete theoretical descriptions have been 
published, for example, references 2 through 4. Thus, in an at 
tempt to obtain a theoretical basis for Owen and Klanfer’s cor 


relation for laminar bubbles, a rather idealized model is postu 


lated, (see Fig. 1 
It is assumed that the flow over this stepped flat plate is in 
compressible and the usual laminar boundary-layer assumptions 


re applicable. It is further assumed that for a small step height, 
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Fic 1. Stepped plate model 





JOURNAL OF THE AERONAUTICAL SCIENCES 


TABLE | 
Laminar Bubble Lengths //6,* 
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the air in the bubble is stagnant and the flow downstream of the 
step and previous to attachment will be essentially a boundary 
layer type flow. The growth of this viscous region is in the direc 
tion of the wall in the manner of a spreading jet. The free-stream 
velocity is assumed to be constant along the entire plate length. 
With this model, the growth of the displacement thickness can 


be predicted on the basis of Blasius’ results, Eq. (1). 


6* = 1.732V/ox/l (1) 
Employing the definition of Rs* and Eq. (1) there results: 
x;/6 = (1/3)Rz.*, %Xa/d.* = (1/3)Ra_* 


Furthermore, since / = x, — x, and with the above, one easily 
»btains Eq. (2) 

1/5,* = (1/3)R5,*[(6a*/5s*)? — 1 (2) 
Originally it was assumed that the separated viscous region grows 
in the wall direction, thus 6,* = 6,* + d and letting d = né6,* 
then 6.* = (nm + 1)6,*. Upon substitution of this relation into 
Eq. (2) the final equation results. 


1/5,* = Rz.*{[(n + 1)? — 1]/3} (3) 


Values of //6,* for several R5,* and a range of »’s are listed in 
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Table 1 with the experimental results from reference 1, for com 
parison 

It is seen from these results that the bubble length grows with 
an increasing step height as did the bubble length in the expe 
mental results for an increasing angle of attack Direct com 
parisons cannot be made since in the experiments diffusiog 
existed in the free stream which resulted in comparatively shorte 
bubbles. A more direct comparison might be made, however, 
if the analysis were based on the growth of the displacement 
thickness over a circular cylinder, for example 

It is also seen that the bubble length grows for an increasing 
R5,* in contrast to the decreasing lengths exhibited in reference] 
for the turbulent type bubbles. This result appears physically 
reasonable if one considers the flow over two plates of different 
initial length at the same free-stream velocity such that the 
separation displacement thicknesses are equal and only the 
viscosities are different, for example, due to a different temper 
ature. In this case, a smaller kinematic viscosity, and thusq 
larger R5,*, would tend to lengthen the separated region due tog 
decrease in the viscous growth in the wall direction. The turby 
lent bubble decreases with increasing R5,* due to the increase ig 
turbulent diffusion of the viscous region. 

While admittedly, this model is a simplification of the actual 
case, certainly such flows exist, and thus it would be interesting 
to test the validity of these results by suitable experiments, 
If reasonable agreement were ascertained, then it is possible that 
this type of boundary-layer analysis could be employed as the 
basis for attacking a more realistic but more difficult case. For 
example, one could employ a relation for the displacement thick- 
ness growth over a curved surface that approximates the free 


stream velocity variation over a separated region 
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